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Abstract 


Consider  a  simulation  experiment  consisting  of  v  independent  vector  observa¬ 
tions  or  replications  across  k  systems,  where  in  any  given  replication  one  and  only 
one  system  is  selected  as  the  best  performer  (i.e.,  it  wins)  based  on  some  perfor¬ 
mance  measure.  Each  system  has  an  unknown  constant  probability  of  winning 
in  any  replication  and  the  numbers  of  wins  for  the  individual  systems  follow  a 
multinomial  distribution.  The  classical  multinomial  selection  procedure  of  Bech- 
hofer,  Elmaghraby,  and  Morse  (Procedure  BEM),  prescribes  a  minimum  number 
of  replications,  denoted  as  v*,  so  that  the  probability  of  correctly  selecting  the  true 
best  system  meets  or  exceeds  a  prespecified  probability.  Assuming  that  larger  is 
better,  Procedure  BEM  selects  as  best  the  system  having  the  largest  value  of  the 
performance  measure  in  more  replications  than  any  other  system. 

In  this  research,  we  use  these  same  v*  replications  across  k  systems  to  form  ( v*)k 
pseudo- replications  that  contain  one  observation  from  each  system,  and  develop 
Procedure  AVC  (All  Vector  Comparisons)  to  achieve  a  higher  probability  of  correct 
selection  (PCS)  than  with  Procedure  BEM.  For  specific  small-sample  cases  and  via 
a  large-sample  approximation  we  show  that  the  PCS  with  Procedure  AVC  exceeds 


n 


the  PCS  with  Procedure  BEM.  In  a  similar  fashion,  we  show  that  with  Procedure 
AVC  we  achieve  a  given  PCS  with  a  smaller  v  than  the  v*  required  with  Procedure 
BEM. 

We  also  consider  the  closely  related  problem  of  estimating  how  likely  each 
system  is  to  be  the  best  under  both  procedures.  Surprisingly,  estimating  pm  is  a 
different  problem  than  estimating  pj,  j  —  1, ...  ,k.  We  show  that  the  variance  of 
the  AVC  estimator  is  never  larger  than  the  variance  of  the  BEM  estimator  (the 
standard  Maximum  Likelihood  Estimator)  and  quantify  the  reduction  in  variance 
with  the  AVC  estimator  for  specific  small-sample  cases  and  asymptotically. 
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CHAPTER  1 


INTRODUCTION 


We  consider  the  problem  of  selecting  the  best  of  a  set  of  alternatives  or  systems 
based  upon  a  comparison  among  them.  To  begin,  we  consider  an  experiment 
consisting  of  a  series  of  independent  trials  across  all  systems  where,  on  any  given 
trial,  one  and  only  one  system  is  selected  as  the  best  performer  (i.e.,  it  wins).  Each 
system  has  an  unknown  constant  probability  of  winning  on  any  trial.  We  count  the 
number  of  wins  for  each  system  over  all  the  trials  and  select  the  system  with  the 
most  wins  as  the  best  system.  In  such  an  experiment,  the  numbers  of  wins  for  the 
individual  systems  follow  a  multinomial  distribution.  The  problem  of  determining 
which  of  the  systems  has  the  largest  probability  of  being  the  best  is  known  as  the 
multinomial  selection  problem  (MSP), 

In  our  context,  we  require  a  quantitative  measure  of  the  performance  of  each 
system  on  each  trial.  Specifically,  we  consider  the  problem  of  determining  which  of 
k  simulated  systems  is  most  likely  to  be  the  best  performer.  A  standard  experiment 
is  to  generate  v  independent  vector  observations  or  replications  (i.e.,  trials)  across 
the  k  systems.  Each  vector  replication  produces  a  single  performance  measure  for 
each  system. 
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Consider  the  following  example.  Suppose  we  are  tactical  war  planning  analysts 
who  are  directed  to  provide  the  Joint  Task  Force  Commander  with  the  best  plan 
to  cripple  the  enemy’s  command  and  control.  “Best”  means  achieving  the  highest 
level  of  cumulative  damage  expectancy  (CDE)  against  a  selected  set  of  targets 
given  current  intelligence  estimates  of  enemy  defense  capabilities  and  available 
friendly  forces.  Our  team  prepares  four  distinct  attack  plans  and  we  simulate  v 
independent  replications  across  all  four  plans.  For  each  replication  we  compare  the 
CDE  among  each  of  the  four  plans.  Since  the  chosen  plan  can  only  be  executed 
a  single  time,  we  select  as  the  best  plan  the  one  that  has  the  largest  CDE  in  the 
most  replications. 

Our  goal  in  an  MSP  is  to  find  the  system  that  is  most  likely  to  be  the  best 
performer  among  the  systems,  as  opposed  to  identifying  the  best  average  performer 
in  the  long  run,  with  a  minimum  amount  of  data.  A  classical  solution  procedure  for 
the  MSP,  Procedure  BEM  (Bechhofer,  Elmaghraby,  and  Morse  1959),  prescribes 
a  minimum  number  of  independent  vector  replications,  v*,  across  all  systems  so 
that  the  probability  of  correctly  selecting  the  true  best  system  meets  or  exceeds  a 
prespecified  probability.  Assuming  that  larger  is  better,  BEM  selects  as  best  the 
system  having  the  largest  value  of  the  performance  measure  in  more  replications 
than  any  other. 

We  propose  a  new  solution  procedure  for  the  MSP  that  requires  no  additional 
data,  but  is  designed  to  increase  the  probability  that  the  best  system  is  so  identi¬ 
fied.  Specifically,  we  propose  using  the  same  v*  replications  across  k  systems  for 
BEM  to  form  ( v*)k  pseudo-replications,  each  containing  one  observation  from  each 
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system,  and  selecting  as  best  the  system  having  the  largest  value  of  the  perfor¬ 
mance  measure  in  more  pseudo-replications  than  any  other  system.  For  specific 
small-sample  cases  and  via  a  large-sample  approximation,  we  show  that  this  new 
procedure,  Procedure  AVC  (All  Vector  Comparisons),  dominates  BEM  in  the  sense 
that  AVC  never  requires  more  independent  replications  than  BEM  to  meet  a  pre¬ 
specified  probability  of  correct  selection  (PCS).  AVC  represents  a  more  efficient 
use  of  the  available  data.  From  a  simulation  design  point  of  view,  we  also  show 
that  by  using  AVC  we  can  achieve  a  given  PCS  with  fewer  replications  than  are 
required  to  reach  the  same  PCS  with  BEM. 

MSP  applications  where  AVC  can  be  applied  include  selecting  the  best  of  a  set 
of  tactical  or  strategic  military  actions  as  presented  earlier.  An  example  in  the 
area  of  structural  engineering  is  finding  the  design  that  performs  best  in  a  one¬ 
time  catastrophic  event,  such  as  an  earthquake.  Some  specific  simulation  examples 
include  selecting  the  schedule  most  likely  to  result  in  completing  all  jobs  on  time; 
selecting  the  investment  portfolio  most  likely  to  provide  the  largest  return;  or 
selecting  the  computer  system  with  the  highest  probability  of  completion  of  a 
series  of  tasks  without  failure.  Each  of  these  applications  involves  the  comparison 
of  quantitative  measures  of  performance  among  competing  systems,  as  opposed  to 
comparing  qualitative  differences.  AVC  requires  a  quantitative  measure  of  system 
performance  for  each  system  in  each  trial  to  be  compared  with  the  performance 
measure  of  other  systems  across  any  or  all  of  the  remaining  trials. 

In  some  MSPs,  the  performance  of  the  alternatives  is  measured  qualitatively, 
rather  than  quantitatively.  Thus,  rather  than  a  comparison  based  on  a  numerical 
measure  of  performance,  we  may  only  know  which  alternative  won  on  any  given 
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trial.  For  example,  in  the  areas  of  marketing  research  or  opinion  surveys,  we 
might  be  interested  in  determining  the  most  popular  brand,  flavor,  etc.,  or  the 
most  favored  candidate  or  position  on  a  political  issue.  In  such  examples,  we 
cannot  compare  the  preference  for  an  individual  alternative  in  one  trial  with  the 
preference  for  one  of  the  remaining  alternatives  in  another  trial.  For  Procedure 
AVC  we  only  consider  quantitative  comparisons. 

In  an  MSP  we  measure  success  in  terms  of  the  probability  we  select  the  true 
best  system,  our  PCS.  However,  the  PCS  contains  no  information  about  how  likely 
the  system  we  selected  as  best  is  to  actually  be  the  best.  In  order  to  obtain  this 
additional  information,  we  can  obtain  a  point  estimate  for  the  probability  of  a 
particular  system  winning  on  any  given  replication.  These  individual  probabilities 
are  multinomial  success  probabilities,  and  we  obtain  different  estimators  for  these 
probabilities  using  BEM  and  AVC.  The  expected  values  of  these  estimators  are 
the  same  for  both  methods;  however,  we  show  that  we  obtain  a  variance  reduction 
with  AVC  over  BEM  and  quantify  this  reduction  for  specific  cases.  In  this  light, 
we  also  compare  AVC  and  BEM  in  terms  of  the  precision  of  the  point  estimators 
provided  by  each  method. 

This  paper  is  organized  into  two  chapters,  each  written  in  the  form  of  a  self- 
contained  journal  article.  Chapter  2  covers  the  multinomial  selection  problem  and 
our  new  approach  to  solving  it,  Procedure  AVC.  Chapter  3  considers  the  related 
point  estimation  problem  and  extends  our  results  into  this  context.  Appendix  A 
contains  proofs,  Appendix  B  presents  additional  simulation  results,  and  Appendix 
C  provides  the  source  code  used. 
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CHAPTER  2 


EFFICIENT  MULTINOMIAL  SELECTION  IN 

SIMULATION 


2.1  Introduction 

Suppose  we  have  k  >  2  independent  populations,  denoted  7Ti ,  7r2 , . . .  ,7r*.  In  a 
simulation  context  each  population  is  a  simulated  system.  We  consider  the  problem 
of  selecting  the  best  of  the  k  systems  based  on  simulated  results  for  all  of  the 
systems. 

Let  Xji  represent  the  ith  replication  from  system  j  of  some  performance  mea¬ 
sure.  Each  system  (7 Tj,  j  =  1, . . . ,  k)  has  an  unknown  constant  probability  (pj,  j  = 
1  of  having  the  largest  value  of  the  performance  measure.  Assume  we  have 

generated  v  independent  replications  from  each  of  the  k  systems.  We  define  the 
best  system  as  the  system  most  likely  to  have  the  largest  performance  measure  in 
any  comparison  across  all  systems.  Such  a  comparison  corresponds  to  a  multino¬ 
mial  trial,  where  one  and  only  one  system  can  win  in  any  given  trial.  Our  objective 
is  to  find  the  system  that  is  most  likely  to  be  the  best  performer  in  a  single  trial 
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among  the  systems,  as  opposed  to  identifying  the  best  average  performer  over  the 
long  run,  with  a  minimum  amount  of  data.  This  is  known  as  the  multinomial 
selection  problem  (MSP). 

A  classical  solution  procedure  for  the  MSP,  Procedure  BEM  (Bechhofer,  El- 
maghraby,  and  Morse  1959),  prescribes  a  minimum  number,  v*,  of  independent 
vector  replications  across  all  systems  so  that  the  probability  of  correctly  selecting 
the  true  best  system  meets  or  exceeds  a  prespecified  probability.  Assuming  that 
larger  is  better,  BEM  selects  as  best  the  system  having  the  largest  value  of  the 
performance  measure  in  more  replications  than  any  other. 

MSP  applications  include  selecting  the  best  of  a  set  of  tactical  or  strategic 
military  actions.  An  example  in  the  area  of  structural  engineering  is  finding  the 
design  that  performs  best  in  a  one-time  catastrophic  event,  such  as  an  earthquake. 
Simulation  examples  include  selecting  the  schedule  most  likely  to  result  in  com¬ 
pleting  all  jobs  on  time;  selecting  the  investment  portfolio  most  likely  to  provide 
the  largest  return;  or  selecting  the  computer  system  with  the  highest  probability 
of  completing  a  series  of  tasks  without  failure.  Each  of  these  applications  involves 
the  comparison  of  quantitative  measures  of  performance  among  competing  systems 
as  opposed  to  comparing  qualitative  measures.  For  the  type  of  MSP  considered 
in  this  study,  we  require  a  quantitative  measure  of  system  performance  for  each 
system  in  each  trial  to  be  compared  with  the  performance  of  other  systems  across 
any  or  all  of  the  remaining  trials. 

Let  X;  =  (AXj,  A2i, . . . ,  A/jj)  represent  the  %  replication  across  all  k  systems. 
Let  Yji  =  1  if  Xji  >  Xu,  for  £  =  1, 2, . . . ,  k,  but  £  ^  j;  and  let  =  0  otherwise. 
In  other  words,  Yjt  —  1  if  X,J%  is  the  largest  observation  in  Xz.  In  case  of  a 
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tie  for  the  largest  value,  we  randomly  select  one  of  the  tied  populations  as  the 
best.  In  terms  of  our  simulation  example  involving  the  best  investment  portfolio, 
Xu  represents  the  return  generated  for  Portfolio  1  on  the  ith  replication.  If  the 
return  on  Portfolio  1  exceeds  the  returns  of  all  the  remaining  portfolios  in  that 
replication,  then  Portifolio  1  wins  (Yu  =  1)  and  all  the  remaining  portfolios  lose 
(Yji  =  0;  j  =  2,3,..  .,&).  If  more  than  one  portfolio  ties  for  the  largest  return 
then  we  randomly  select  one  of  these  to  be  the  winner  in  that  replication. 

Suppose  that  there  are  v  independent  replications  across  all  systems,  and  let 
Yj  =  J2i=i  Yji  represent  the  number  of  times  system  j  wins  out  of  these  v  repli¬ 
cations.  Let  Pj  =  Pr{Arjj  >  Xu,  Ml  ^  j}  where  0  <  pj  <  1  and  J2j=iPj  =  L 
Then  X^=1  Yj  —  v  and  the  k-v ariate  discrete  random  variable  Y  =  (Y,  Y2, . . . ,  Yk) 
follows  a  multinomial  distribution  with  success  probabilities  p  =  (pi,p2,  ■  ■  ■  ,pk)- 
Therefore,  the  probability  mass  function  for  Y  with  parameters  v  and  p  is 


vl 


Pr{Y  =  j/i,  Y2  =  y2,...,  Yk  =  yk)  =  j  II  Pj  • 


UKj=1yj\ 


3= 1 


For  a  standard  MSP  with  v  replications  from  each  system,  we  have  v  indepen¬ 
dent  comparisons  (trials)  to  use  in  selecting  the  best  system.  For  an  experiment 
involving  physical  measurements  of  system  performance  at  common  times,  it  makes 
sense  to  group  the  observations  across  systems  for  the  same  trial  due  to  possible 
variations  in  extraneous  factors  beyond  the  experimentor’s  control.  Also,  in  some 
instances,  the  performance  of  the  competing  systems  may  be  measured  qualita¬ 
tively,  or  may  only  indicate  which  system  won  in  a  given  trial.  In  such  cases,  there 
is  no  quantitative  measure  we  can  compare  across  systems  in  different  replications. 
Examples  include  marketing  research  or  opinion  surveys  where  the  data  collected 
are  categorical  (e.g.,  favorite  brand,  food,  or  political  candidate). 
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Due  to  convention  and  convenience  when  comparing  simulated  system  re¬ 
sponses,  the  responses  are  typically  grouped  by  replication,  corresponding  to  a 
trial  in  a  physical  experiment.  Grouping  system  responses  in  this  fashion  is  ar¬ 
bitrary  and  since  our  simulated  responses  are  quantitative,  we  can  compare  any 
observation  from  one  system  with  any  observation  from  each  of  the  remaining  sys¬ 
tems.  This  means  that  a  single  observation  from  system  1  can  be  grouped  in  a 
vector  comparison  with  any  one  of  the  v  observations  from  system  2,  and  with  any 
one  of  the  v  observations  from  system  3,  and  so  on,  up  to  and  including  any  one 
of  the  v  observations  from  system  k.  Since  there  are  v  observations  from  system  1 
as  well,  this  gives  us  a  total  of  vk  vector  comparisons  (trials)  that  can  be  formed 
with  v  independent  observations  from  the  k  systems.  We  incorporate  this  setup 
in  a  new  MSP  procedure,  which  we  call  AVC,  for  All  Vector  Comparisons.  By 
performing  only  the  v  vector  comparisons  where  the  observations  for  each  system 
are  from  the  same  replication,  as  is  done  with  BEM,  we  disregard  the  information 
available  from  the  remaining  vk  —  v  comparisons. 

Our  results  suggest  a  number  of  advantages  of  AVC  over  BEM.  For  specific 
small-sample  examples,  we  show  that  AVC  has  a  larger  probability  of  correct  se¬ 
lection  (PCS)  than  BEM  for  a  fixed  v.  We  show  this  analytically  for  small  values 
of  v  and  k,  and  also  present  simulation  results  for  up  to  k  =  10  systems  and  v  =  50 
vector  replications.  Looking  at  these  results  from  a  slightly  different  perspective, 
we  also  demonstrate  acheivement  of  a  desired  PCS  with  a  smaller  value  of  v  when 
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using  AVC  as  compared  to  BEM.  The  first  perspective  emphasizes  a  more  effi¬ 
cient  use  of  the  available  data  to  increase  PCS.  The  second  view  points  towards  a 
more  efficient  way  to  design  a  simulation  experiment  using  the  smallest  value  of  v 
required  to  achieve  a  desired  PCS. 

Unlike  BEM,  the  PCS  for  AVC  depends  on  the  distributions  of  the  simulation 
outputs,  not  just  on  pi, . . .  ,pk-  However,  we  also  show  that  the  dependence  is  weak. 
This  fact,  along  with  the  difficulty  of  analytically  evaluating  the  PCS  of  AVC  for 
even  small  k  and  v,  leads  us  to  a  large-sample  approximation  (LSA)  for  the  PCS 
using  AVC.  As  v  — >■  oo,  any  distributional  differences  in  PCS  with  AVC  disappear. 
Therefore,  our  LSA  is  distribution  independent  and  converts  an  AVC  problem  into 
an  equivalent  BEM  problem.  Our  LSA  demonstrates  that  asymptotically  the  PCS 
with  AVC  is  larger  than  the  PCS  with  BEM.  Additionally,  this  LSA  shows  that 
AVC  can  provide  better  discrimination  between  the  systems  at  the  same  level  of 
confidence  and  with  the  same  data. 

This  paper  is  organized  as  follows:  We  first  provide  a  brief  review  of  the  MSP 
and  the  classical  approach  to  solving  it.  Then  we  describe  our  new  procedure,  AVC, 
and  present  analytical  results  covering  a  variety  of  specific  population  distributions 
for  the  performance  measures.  Our  LSA  is  then  presented  by  recasting  PCS  in 
terms  of  a  point  estimation  problem  for  the  multinomial  success  probabilities, 
Pj,  j  =  1  Empirical  results  follow  for  specific  distributions  and  include 

simulations  designed  to  test  the  robustness  of  our  LSA. 
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2.2  Background 


Bechhofer,  Elmaghraby  and  Morse  (1959)  describe  a  single-stage  procedure  for  se¬ 
lecting  the  multinomial  event  (population  or  system)  which  has  the  largest  success 
probability.  BEM  requires  the  specification  of  P*  (where  l/k  <  P*  <  1),  a  min¬ 
imum  probability  of  correctly  identifying  the  population  with  the  largest  success 
probability  (i.e. ,  the  best  population),  and  9*  (where  1  <  9*  <  oo),  the  minimum 
ratio  of  the  largest  success  probability  to  the  second  largest  success  probability  that 
we  want  to  be  able  to  detect.  The  procedure,  as  adapted  to  simulation,  consists 
of  the  following  steps: 

Procedure  2.1  (BEM) 

1.  For  given  k  and  9* ,  find  the  minimum  value  ofv,  denoted  v* ,  that  guarantees 
that  the  PCS  is  at  least  P* . 

2.  Generate  v*  independent  replications  for  each  population. 

3.  Compute  Yj  =  Yfl-i  Yji,  for  j  —  1, 2, . . . ,  k. 

4 .  Let  Y(1)  <  F( 2)  <  •  •  •  <  Y(k)  be  the  ranked  sample  counts  from  step  3.  Select 
the  population  associated  with  the  largest  count,  Y^,  as  the  best  population. 
In  case  of  a  tie  for  the  largest  count,  randomly  select  one  of  the  tied  popula¬ 
tions  as  the  best. 

To  determine  the  appropriate  v*  in  step  1,  let  ppj  <  p[2]  <  •  •  •  <  p[k\  denote  the 
ranked  success  probabilities  for  the  k  populations.  Since  only  values  of  the  ratio 
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9  =  p[k]/p[k-i]  greater  than  or  equal  to  9*  are  of  interest,  we  are  indifferent  between 
the  best  and  the  next-best  population  for  values  of  9  <  9*.  A  procedure  of  this 
type  is  referred  to  as  an  indifference-zone  approach.  Select  v*  as  the  minimum 
number  of  independent  vector  observations  required  to  achieve  a  PCS  greater  than 
or  equal  to  P*  whenever  9  >  9*. 

We  define  the  least  favorable  configuration  (LFC)  of  p  =  (ppppp],  •  •  •  ,P[k])  as 
the  configuration  where  PCS  is  a  minimum  over  all  configurations  with  9  >  9* 
(Gibbons,  Olkin,  and  Sobel  1977).  If  we  obtain  a  PCS  >  P*  with  our  selected  v* 
under  the  LFC,  then  a  PCS  of  at  least  P*  can  be  guaranteed  for  any  configuration 
of  p  with  9  >  9*.  Keston  and  Morse  (1959)  prove  that  the  LFC  for  BEM  is  given 
by 

l 

Pm  =  Pm  =  ■■■  =  P[k- 1]  =  0*  +  k  _  1 

Q * 

P[k]  =  9*  +  k  -  1 '  ^ 


Although  we  only  need  to  consider  the  LFC  for  designing  sampling  plans,  the  PCS 
can  be  calculated  for  any  p  with  p^}  >  P\k-\)  as  follows. 

Let  T\[j]  be  the  population  associated  with  p^  and  let  represent  the  number 
of  wins  for  7ryj.  Thus,  the  subscripts  for  the  populations  and  the  associated  number 
of  wins  are  based  on  the  ranking  of  the  pj s.  We  refer  to  the  PCS  using  BEM  for 
a  fixed  k  and  v  as  PCSbera.  For  any  fixed  k  and  v,  PCSbemcan  be  expressed  as 


PCSbem(p)  =  £ 
y 


1  w! 

W)  n‘„  m  ■ 


where  the  summation  is  over  all  vectors  y  =  (?/[i],  •  •  • ,  ?/[fe])  such  that  J2j=iUj  = 
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V,  V[k]  >  y\j]  {j  =  1,2  1),  where  t(  y)  is  a  function  of  y[x  j, . . . ,  y[k],  repre¬ 

senting  the  number  of  populations  tied  for  the  most  wins  (Bechhofer,  Elmaghraby, 
and  Morse  1959). 


2.3  All  Vector  Comparisons  (AVC) 


We  propose  a  method  to  provide  a  PCS  greater  than  or  equal  to  PCSbem  (in  at 
least  some  cases)  using  the  same  replications  X,,  i  =  1, 2, . . . ,  v.  We  use  the  BEM 
parameters  k ,  P* ,  and  6* ,  and  we  execute  the  first  step  of  BEM  to  find  a  value 
of  v* .  However,  rather  than  comparing  the  ith  replication  for  each  system  with 
the  ith  replications  of  the  other  systems,  consider  instead  a  total  of  ( v*)k  pseudo¬ 
replications  formed  by  associating  each  X}i  (j  —  1,2,  i  =  1,2,  ...,v*), 

with  all  possible  combinations  of  the  remaining  X^h  (t  =  1, 2, . . . ,  k\  £  ^  j;  h  = 
1,2, ...  ,v*).  Each  such  pseudo-replication  contains  one  observation  from  each 
population.  Note  that  the  ( v*)k  pseudo- replications  include  the  v*  independent 
replications  from  which  the  pseudo-replications  are  formed. 

Define 

n  <Kxn  -  Xu,)  (2-2) 

aj  =1  rip  =  1  dk~  1  £=l]£^j 

for  j  =  1,  2, . . . ,  A:  with 


</>(a) 


'  1,  a  >  0 

0,  a  <  0 

randomly  assign 
0  or  1,  a  =  0. 


Thus,  Zj  represents  the  number  of  times  out  of  vk  pseudo-replications  that  pop¬ 
ulation  7 Tj  wins  (ties  broken  randomly)  and  Xq=i  Zj  =  vk. 
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As  a  specific  illustration  of  how  the  pseudo-replications  are  formed,  consider 
k  =  3  systems  with  v  —  2  observations  generated  for  each.  Our  original  replications 


and  counts  (Yj,  j  =  1, 2, 3)  are  then 

Xu  X21  X3l 
X 12  X‘22  X32 
~Yi  %  YT 

We  will  have  a  total  of  vk  (23  =  8)  pseudo-replications  including  our  two  original 
replications  when  using  AVC.  These  pseudo- replications  and  the  associated  counts 
{Zj,  j  =  1,2,3)  are 


AT 

A2i 

A3i 

An 

X2\ 

a32 

An 

X22 

A3i 

An 

A22 

a32 

Xn 

X2\ 

A3i 

Xn 

X2\ 

a32 

AT2 

A22 

A3i 

at2 

a22 

a32 

Z  i  Z2  Z3 


Our  new  procedure  consists  of  the  following  steps: 


Procedure  2.2  (AVC) 

1.  Given  values  for  k,  P* ,  and  9* ,  use  step  1  of  Procedure  BEM  to  determine  a 
value  for  v* . 

2.  Generate  v*  independent  replications  for  each  population  and  construct  the 
additional  (v*)k  —  v*  pseudo-replica, tions  possible  with  one  value  from  each  of 
the  populations. 

3.  Compute  Zj  using  Equation  (2.2). 
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4 .  Let  Z(i)  <  Z( 2)  <  •  •  •  <  Zfk)  be  the  ranked  sample  counts  from  step  3.  Select 
the  population  associated  with  the  largest  count,  Z^,  as  the  best  population. 
In  case  of  a  tie  for  the  largest  count,  randomly  select  one  of  the  tied  popula¬ 
tions  as  the  best. 

Later  we  demonstrate  for  specific  cases  that  the  PCS  with  AVC,  referred  to  as 
PCSavc,  is  greater  than  or  equal  to  PCSbem. 

As  written,  step  1  of  Procedure  AVC  uses  the  same  number  of  replications  as 
BEAL  Suppose  we  modify  step  1  to  use  the  minimum  v  where  PCSavc  >  P*.  We 
demonstrate  later  that  a  smaller  number  of  replications  are  required  with  AVC 
relative  to  BEM  to  achieve  P*.  We  provide  such  values  of  v  in  this  paper. 

PCSavc  can  be  expressed  as 

PCSavc(p)  =  £^Pr{Zw  =  ^|,...,Zw  =*!*]}, 

(z) 

where  the  summation  is  over  all  vectors  z  =  (^q, . . . ,  Z[k])  such  that  J2j=i  Zj  — 
vk,  Z[k\  >  %],  j  =  1,  2, . . . ,  k  —  1,  where  t( z)  is  a  function  of  Z[i], . . . ,  z^],  repre¬ 
senting  the  number  of  populations  tied  for  the  most  wins.  Unfortunately,  Z  does 
not  follow  a  multinomial  distribution,  so  we  must  refer  to  the  distributions  of  the 
original  observations,  XJly  to  calculate  PCSavc.  Analytical  and  simulation  results 
using  a  number  of  different  population  distributions  show  that  PCSavc  depends 
weakly  on  the  underlying  distributions  of  the  XjiS. 


14 


2.4  Analytical  Results 


The  following  analytical  study  illustrates  a  number  of  important  properties  of  the 
AVC  method.  First,  we  demonstrate  the  improvement  possible  with  AVC  for  spe¬ 
cific  cases.  We  also  show  a  weak  dependence  in  the  AVC  results  on  the  underlying 
population  distributions  for  the  Xji.  Lastly,  we  demonstrate  the  difficulty  in  ob¬ 
taining  analytical  results  for  even  a  small  number  of  populations  and  observations, 
and  thus  provide  motivation  for  our  large  sample  approximation  of  PCSavc  which 
is  distribution  independent. 

2.4.1  Small-Sample  Results 

Initially,  we  restrict  our  attention  to  continuous  distributions  for  the  XjiS,  which 
eliminates  the  possibility  of  ties  among  the  observations.  We  let  7^  be  the  best 
population  and  assume  all  the  remaining  populations,  7T[i], . . . ,  7r[*_i],  are  identically 
distributed.  This  setup  gives  us  the  LFC  for  BEM.  We  also  consider  all  population 
distributions  to  belong  to  the  same  parametric  family.  We  will  calculate  PCSavc  by 
conditioning  on  the  joint  density  of  all  the  order  statistics  for  the  v  independent 
replications  from  7^]. 

Consider  a  set  of  v  vector  replications  across  all  populations.  Combine  all  the 
observations  from  all  populations  and  rank  them  from  smallest  to  largest.  Refer 
to  each  observation  by  its  rank  and  consider  permutations  of  these  ranks.  For  any 
such  permutation  we  can  determine  the  value  of  and  calculate  the  probability 
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of  obtaining  that  arrangement  of  ranks.  We  refer  to  such  an  arrangement  as  a  rank 
order.  Recall  that  represents  the  number  of  times  the  best  population,  7 r[fe], 
wins  out  of  the  vk  pseudo-replications.  For  illustrative  purposes,  let  X  represent  an 
observation  from  7T[fc]  and  let  O  represent  an  observation  from  any  of  the  remaining 
inferior  populations. 

As  an  example,  suppose  k  =  3,  v  =  2.  Then 

Pr{^[3]  =  8}  =  Pr{0(p  <  0(2)  <  0(3)  <  0(4)  <  A"(p  <  X(2)}  (2.3) 

Pr{%]  =  6}  =  4  Pr{0(i)  <  0(2)  <  0(3)  <  A"(i)  <  0(4)  <  A"(2)}  (2.4) 

Since  we  do  not  know  which  of  the  four  Os  is  associated  with  which  rank,  we 
must  account  for  all  permutations  of  the  Os  that  result  in  a  different  combination 
of  adjacent  Os.  For  probability  statement  (2.3),  there  is  only  one  combination  of 
adjacent  Os  from  the  4!  permutations  of  the  Os  that  is  less  than  both  Xs.  In 
the  rank  order  for  probability  statement  (2.4),  since  any  one  of  the  Os  can  be 
associated  with  0(4),  we  have  four  distinct  combinations  (in  terms  of  which  set  of 
Os  are  adjacent)  that  result  in  this  one  rank  order.  This  is  why  the  coefficient 
‘4’  appears  on  the  right-hand  side  of  Equation  (2.4).  In  general  this  coefficient  is 
,  where  n  =  v(k  —  1)  is  the  total  number  of  observations  from  the  inferior 
populations  and  r  is  the  largest  number  of  these  observations  that  are  adjacent. 
Similar  arguments  can  be  used  to  derive  expressions  for  possible  values  of  Z\k]  for 
integers  k,  v  >  2.  For  this  example,  there  is  only  a  single  rank  order  that  results 
in  each  value  of  Z[ky  As  k  or  v  get  even  moderately  large,  there  will  be  multiple 
rank  orders  that  result  in  the  same  value  for  Z^y  In  addition,  for  larger  k  and 
v  we  must  also  take  into  consideration  how  many  of  the  combinations  contain  an 
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observation  from  each  of  the  inferior  populations.  Therefore,  the  calculation  of  the 
probability  of  each  value  of  Z[k]  becomes  extremely  tedious  with  increasing  k  or  v. 

Restricting  our  attention  to  k  —  2  populations,  it  is  interesting  to  note  that 
the  vector  comparisons  with  AVC  are  analagous  to  the  comparisons  that  form  the 
Wilcoxon  rank-sum  statistic  (Randles  and  Wolfe  1979).  Let  W  equal  the  sum  of 
the  ranks  of  the  observations  from  the  best  population.  Then  W  is  the  Wilcoxon 
rank-sum  statistic  and  our  Z[2]  is  the  Mann- Whitney  U  statistic.  Therefore,  W 
can  be  expressed  as  a  function  of  our  Z[2]  as 

W  =  zm+l-(v  +  l). 

In  terms  of  W,  AVC  always  makes  a  correct  selection  for  W  >  E[W]  (incorrect 
selection  for  W  <  E[W]),  where  E[kF]  is  the  expected  value  of  W  under  the 
assumption  that  the  two  populations  are  identical  in  distribution. 

If  we  specify  a  particular  distribution  family  for  our  populations,  then  we  can 
derive  formulas  to  compare  PCSavc  with  PCSbem  for  small  k  and  v.  We  present 
results  for  exponential,  continuous  uniform,  and  Bernoulli  distributions. 
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2.4.2  Exponential 


First,  consider  X  ~  exp(A)  and  O  ~  exp(/r)  and  let  A  <  //,,  where  A  >  0  and 
/i  >  0  are  exponential  rates.  This  particular  example  was  suggested  by  Goldsman 
(1995)  in  some  early  discussions  concerning  this  research.  For  k  ==  2,  v  =  2,  we 
have  p[2]  —  Pr{A'  >  0}  =  ///(A  +  //)  and  pp]  =  Pr{A"  <  0}  =  A/(A  +  //).  To 
calculate  PCSbem,  we  need  to  consider  vectors  Y  =  (Ip],  1[2])  such  that  Fp]  >  Fpp 
With  v  =  2,  the  only  possible  Ys  with  Fp]  >  Ip]  are  (0,2)  and  (1, 1).  This  gives 
us 

PCSbem  =  Pr{F[2]  =  2}  +  ^  Pr{F[2]  =  1}  (2.5) 

=  P[2]  +  l2P[l}P{2] 

F 

A  +  fi 

Similarly,  to  calculate  PCSavc,  we  need  to  consider  vectors  Z  =  (ZppAp])  such 
that  Z[2]  >  Z[ij.  With  vk  =  4,  the  only  Zs  with  Zp ]  A  Z[i]  are:  (0,4),  (1,3),  (2,2). 
So 

PCSavc  =  Pr{Z[2]  =  4}  +  Pv{Z[2]  =  3}  +  ^  Pr{Z[2]  =  2}.  (2.6) 


When  A”  ~  exp(A),  the  joint  distribution  of  (ATp],  Ar[2])  is  fn(a,  b )  =  2A2e_A^a+6h 
The  probabilities  on  the  right-hand  side  of  Equation  (2.6)  can  then  be  found  as 
follows: 


18 


Pt{Z[2]  =  4}  =  Pr{0{1)  <  0(2)  <  X(1)  <  X(2)} 


nb 

(1  -  e~tia)22X 

. 


2e~xP+b')d,a  db 


(2A  +  /i)(A  +  ji) 


nb 

(1  -  e-^)(e-^  -  e-^)2AVA(a+b)dad& 


2Ap2 

(2A  +  p)(A  +  fi )2 


n6 

(e_/ia  -  e-^bf2\2e-x{a+b)dadb  + 

2 ii^1  “  ^“K^AV^+^da  db 


2A/x(A2  +  4A//  +  p2) 
(2A  +  p)(A  +  2//)  (A  +  p)2 


Therefore, 


PCSavc  = 


p( A2  6Ap  +  2/i2 


(2A  +  p)(A  +  2p)(A  +  p) 

Given  expressions  for  both  PCSbem  and  PCSavc,  we  can  find  the  increase  in 


PCS  with  AVC  as 


APCS  =  PCSavc  -  PCS1 


=  _ A ^  -  A> _ >0  (27) 

(2A  +  //)(A  +  2/,)(A  +  ^)  k  1 

The  (p  —  A)  >0  term  in  Equation  (2.7)  shows  that  when  X  is  the  best  population, 

AVC  always  shows  an  improvement  in  PCS  over  BEM. 

Similar  calculations  for  k  =  2  and  v  =  3  result  in 


bem  P2(  3A  +  //) 


PCS 


PCS avc  = 


(A  +  p)3 

p2(20A5  +  159AV  +  344A3p2  +  273AV  +  92Ap4  +  12p5) 
(3A  +  2p)(2A  +  3p)(2A  +  p)(A  +  2p)(A  +  p)3 
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and 


=  A2//2 (16 A2  +  37  Xp  +  16/i2)(/x  -  A) 

(3A  +  2p)(2\  +  3/i)(2A  +  /r)(A  +  2/r)(A  +  yu)3 

Also  for  k  =  3  and  w  =  2  we  obtain 


PC'S4' 


PCS'8 


2p2 

(A  +  2/i)(A  +  /r) 

2/z3(3A4  +  100A3/r  +  267A  V  +  194V  +  36/P) 
(2A  +  3//)(2A  +  //)( A  +  4/<)(A  +  3//)(A  +  2//)(A  +  /r)2 


APC  9  —  2A/P(4A3  +  41  A2//  +  84A/P  +  41/i3)(//  —  A) 

(2A  +  3/u)(2A  +  /i)(A  +  4/r)(A  +  3/r)(A  +  2/i)(A  +  /i)2 

Each  of  the  expressions  above  for  the  difference  in  PCS  between  AVC  and  BEM 
has  a  (/i  —  A)  >  0  term,  implying  an  increase  in  PCS  with  AVC. 


2.4.3  Continuous  Uniform 


Consider  k  =  2  and  v  =  2  where  Ar  and  O  come  from  continuous  uniform  distri¬ 
butions.  Specifically  X  ~  17(0,  P)  and  O  ~  17(0,  A),  where  0  <  A  <  P.  Then  we 
have  pp]  =  Pr{A  >  0}  =  (2 B  —  A)/2B  and  pqj  =  Pr{A  <  0}  =  A/2B.  From 
Equation  (2.5)  we  obtain 


PCSbem  - 


2B- A 


The  joint  distribution  of  (XppXp])  is  fn(g,h )  =  2/P2.  Then  using  Equa¬ 
tion  (2.6)  we  find 


PCS’avt  = 


A2  +  2AP  -  6 P2 
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2.4.4  Bernoulli 

As  an  illustration  of  how  AVC  compares  to  BEM  for  discrete  distributions,  let 
X  ~  Bern(pz)  and  O  ~  Bern(p0)  with  px  >  p0.  For  k  =  2,  v  =  2,  we  have 
P[ 2]  =  Pr{A"  >  0}  =  (px  +  1  -  p0)/2  and  pw  =  Pr{X  <  0}  =  (p0  +  1  -  px)/2. 
From  Equation  (2.5)  we  obtain 

PCSbern  =  ^(px  +  1  -  p0). 

PCSavc  can  be  calculated  in  a  similar  fashion  as  we  did  for  our  exponential  and 
continuous  uniform  examples, 

PCSaVC  =  ~(pI  -  pi)  +  \pxPo{Px  -  Po)  +  \{Px  -  Po)  + 
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Then  calculating  the  difference  in  PCS, 


APCS  =  ^(1  +  2 pxp0  -px-  p0){Px  ~  Po )  >  0.  (2.8) 

We  again  see  a  term,  (px  —  p0)  >  0  hr  Equation  (2.8),  which  illustrates  an  improve¬ 
ment  in  PCS  with  AVC  when  X  is  the  better  population. 


2.4.5  Effect  of  Population  Distributions  on  AVC 

Our  analytical  results  are  presented  for  single  values  of  k  and  v  for  each  of  the 
distributions.  To  compare  PCSavc  across  the  distributions,  we  use  a  fixed  value 
of  k  and  v  and  vary  9  =  p\k}/p[k- 1]-  We  display  these  results  using  the  ratio  of 
improvement  PCSavc/PCSbem  plotted  against  9.  Recall  that  PGSbem  is  distribution 
independent. 

Figure  2.1  shows  results  for  exponential,  continuous  uniform,  and  Bernoulli 
populations  at  k  =  2,  v  =  2.  The  exponential  and  uniform  curves  depend  only 
on  9  and  approach  1  as  9  approaches  infinity.  However,  for  the  Bernoulli  results, 
we  must  specify  a  value  for  either  px  or  p0,  which  forces  an  upper  limit  on  the 
value  of  9  as  px  approaches  1  or  as  pQ  approaches  0.  We  set  px  =  7/8  and  then 
limPo_>o0  =  15.  Figures  2.2  and  2.3  display  results  for  exponential  and  continuous 
uniform  populations  at  k  =  2,  v  =  3  and  k  —  3,  v  =  2  respectively. 

There  are  some  significant  results  to  gleen  from  Figures  2.1-  2.3.  First,  notice 
the  spread  between  the  results  for  the  different  populations  is  quite  small  (a  maxi¬ 
mum  of  roughly  5%  for  k  =  3,  v  =  2,  and  on  the  order  of  1%  for  k  =  2,  v  —  2  and 
k  =  2,  v  =  3),  and  this  spread  drops  off  as  9  approaches  one.  This  indicates  that 
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1.08 


Figure  2.1:  Exponential,  Uniform  and  Bernoulli  Populations:  k  =  2,  v  =  2,  px  =  | 


as  the  difference  between  the  best  and  the  next  system  gets  smaller,  the  distribu¬ 
tional  dependence  becomes  weaker.  Also  notice  that  the  maximum  improvement 
with  AVC,  indicated  by  the  peaks  in  the  curves,  occurs  roughly  over  a  range  of  9 
between  1.5  and  3  for  all  the  distributions.  This  covers  most  of  the  practical  range 
of  9  included  in  standard  tables  and  used  by  experimentors. 

The  weak  distributional  dependence  of  PCSavc,  along  with  the  difficulty  of 
computing  PCSavc  for  small  k  and  v,  motivates  the  following  large-sample  approx¬ 
imation  (LSA). 
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Figure  2.2:  Exponential  and  Uniform  Populations:  k  —  2,  v  =  3 


2.5  Large- Sample  Approximation 

The  results  presented  so  far  for  small  k  and  v  show  that  PCSavc  is  weakly  distribu¬ 
tion  dependent.  By  redefining  our  PCSbem  and  PCSavt  in  terms  of  point  estimators 
for  each  of  the  individual  system  success  probabilities,  we  arrive  at  distribution-free 
results  as  the  sample  size  goes  to  infinity. 

2.5.1  Preliminaries 

Using  our  previous  notation  we  have 

pj  =  Pr{Xji>Xti,W?j}. 
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Figure  2.3:  Exponential  and  Uniform  Populations:  k  =  3,  v  =  2 


Let  the  distribution  of  Xjt  depend  upon  the  sample  size,  Xj  ~  Fjv\  We  construct 
the  FjV^  such  that  the  F converge  to  a  common  distribution,  F,  for  all  j  as  v 
approaches  infinity  and  for  finite  v 


Pr  {Xji  >  Xje  |  sample  size  v} 


Define 


Pj(v) 


(k-l)6 

V* 

5 

Vi 


j  =  1 
j  +  I- 


(2.9) 


Yj(y)  =  number  of  wins  for  system  j  under  BEM  with  sample  size  v 
Zj(v)  =  number  of  wins  for  system  j  under  AVC  with  sample  size  v 

which  gives  us  point  estimators 
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p3(v)  = 


Zj(v)' 
vk 

So  our  BEM  estimators  are  denoted  by  pj  and  our  AVC  estimators  by  pj.  Under 
(2.9)  Population  1  is  the  best.  Notice  that  (ignoring  the  asymptotically  vanishing 
probability  of  a  tie) 

PCS'”"  =  Pr{ft  >  ft,  Vi  /  1} 

PCS~  =  Pr{ft>ft,Vjftl}. 

Our  approach  is  based  on  the  fact  that  standardized  versions  of  p  and  p  are 
asymptotically  multivariate  normal  (MVN).  However,  when  the  distributions  are 
fixed  then  as  the  sample  size  increases,  both  PCSbem  and  PCSavc  approach  1, 
masking  the  differences  between  the  two  procedures.  To  eliminate  this  effect  and 
isolate  the  improvement  with  AVC,  wfe  simultaneously  let  the  ratio  of  the  differ¬ 
ences  between  the  populations  approach  1  at  the  canonical  rate  of  1/ y/v,  as  shown 
in  (2.9). 

2.5.2  BEM  Estimators 

Consider  the  asymptotic  behavior  of  PCSbem  as  the  number  of  vectors,  v,  goes  to 
infinity.  Our  approach  is  structured  around  a  result  presented  by  Lehmann  (1986, 
pp.  478-479)  which  we  state  below  as  a  lemma. 

Lemma  2.1  (Asymptotic  Distribution  of  Standardized  BEM  Estimators )  Let  Y  = 
(l'i (v),Y2(v),  . . .  ,Yk{v ))  be  distributed  as  a  multinomial  random  variable  with  pa¬ 
rameters  v  and  p  =  (pi(v),p2(v), . . .  ,pk(v)).  Define  Pj(v)  as  in  equation  (2.9). 
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Let 


j  =  1,2  ,...,k. 


Then  as  v  — »  oo 


Ai(v) 


Yj(v)  -  v/k 


(  A2(v)  ^ 

V  Mv) 


'■ 

/  l//c(l  -  1/k) 

-i/k(i/k)  ■ 

■  — i/*(i A)  \ 

•  - 

MVN 

-iA(iA) 

l/k(l  —  1/k)  ■ 

■  -l/k(l/k) 

[s) 

V  -i/k(i/k) 

-l/k(l/k)  ■ 

■  1/*(1  -  1/k)  J 

Proof:  See  Lehmann  (1986). 

Lemma  2.1  is  critical  to  proving  the  following  theorem  which  we  later  use  to 
equate  asymptotic  PCSbem  with  asymptotic  PCSavc. 


Theorem  2.1  (Asymptotic  PCSbem) 

Let  Y  =  (Yi(v),Y2(v),  ... ,  Yk(v))  be  distributed  as  a  multinomial  random  vari¬ 
able  with  parameters  v  and  p  =  (pi{v),p2(v), . . .  ,pk(v)),  with  Pj  defined  as  in 
equation  (2,9).  Then 


bem 


PCS  =  Jhn  pr{lri('e)  >  Tj(n)}  =  Pr 


where 


(Q2 ) 

(o) 

/ 

~  MVN 

V  Qk  ) 

v0 ) 

V 

Proof:  See  Appendix  A. 


1/2  \  ’ 
1  )  . 


(2.10) 
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2.5.3  AVC  Estimators 


Consider  the  asymptotic  behavior  of  PCSavc  as  the  number  of  vectors,  v,  goes  to 
infinity.  Our  approach  is  structured  around  a  result  presented  by  Lehmann  (1963, 
pp.  964-965)  and  also  in  Randles  and  Wolfe  (1979,  p.  107).  We  state  this  result 
below  as  Lemma  2.2. 

We  notice  that  our  AVC  estimator  is  a  k- sample  U-statistic,  where  p  = 
(pi,P2,  ■  ■  ■  ,Pk)  represents  our  parameter  for  the  system  success  probabilities.  From 
Randles  and  Wolfe  (1979,  p.  104),  we  say  p  is  estimable  of  degree  (1, . . . ,  1)  for 
distributions  (F\, . . . ,  Fk)  in  some  family  F,  if  (1, . . . ,  1)  are  the  smallest  sample 
sizes  for  which  there  exists  an  unbiased  estimator  of  p  for  every  (F\, . . . ,  Fk)  €  T . 
Formally  stated  we  have 

E{Fu...,Fk)[h('j\X ii, . . . ,  Xfci)]  =  pj 

for  j  =  1, . . . ,  k  and  we  define  our  ^-sample  symmetric  kernel  h^\-)  as 

k 

hU}  =  n  —  Xu).  (2.11) 

So  we  have  =  1  if  the  observation  from  the  jth  system  is  the  largest  in  any 
vector  comparison  across  all  systems.  Substituting  our  kernel  from  (2.11)  into 
(2.2)  and  dividing  by  vk,  we  have  presented  our  pj ,  j  =  1,  2, . . . ,  k  as  fc-sample 
U-statistics.  We  can  then  state  the  following  lemma,  where  some  of  the  notation 
has  been  simplified  to  match  the  context  of  our  problem. 

Lemma  2.2  (Asymptotic  Distribution  of  AVC  Estimators)  Let  pi,...,pk  be  k- 
sample  (k  >  1 )  U-statistics,  with  pi  corresponding  to  a  parameter  pi  of  degree 
(1, . . . ,  1)  and  symmetric  kernel  for  i  =  1, . . . ,  k.  Let  N  =  kv,  where  v  is 
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the  sample  size  from  each  of  k  populations.  Then  the  joint  limiting  distribution  of 


(  y/N(pi  ~  Pi )  ^ 

(  0  \ 

=>  MVN 

,  E  =  (a(a’ft)) 

\  \fN (Pk  -  Pk)  y 

^0  J 

where 


for  A i  =  limjv^oo  v/N. 


aM)  =  £  1  f  W) 

i=  1 

The  quantities  are  given  by 


d"'*’  =  Cov\h}?,h£ 
=  e[h^hI: 


m 

i2 


PaPb 


where 


(2.12) 


=  /.W(A'w„A%,...,A'iA) 

and  the  sets  (ai,  ct2, . . . ,  a*, . . . ,  a k)  and  (Pi,  P2,  ■  ■  • ,  Pi, . . . ,  Pk)  have  only  the  ith 
element  in  common,  where  the  elements  in  each  set  represent  positive  integers  and 
a,b  G  k}. 


Proof:  See  Lehmann  (1963). 

Proceeding  in  much  the  same  manner  as  we  did  in  moving  from  Lemma  2.1  to 
Theorem  2.1,  we  can  define  PCSavc  as  a  probability  statement  involving  a  function 
of  k  and  the  maximum  of  (k  —  1)  random  variables  created  by  subtracting  one  of 
the  standardized  pjS  from  each  of  the  remaining  standardized  fps,  i  /  j.  These 
random  variables  are  the  Q2,  Q.s,  ■  ■  . ,  Qk  in  the  following  theorem. 
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Theorem  2.2  (Asymptotic  PCSavcJ 


Letpi,...,pk  be  k-sample  (k  >  1)  U-statistics,  with  pi  corresponding  to  a 
parameter  Pi  of  degree  (1, . . . ,  1)  and  symmetric  kernel  for  i  =  1 , ,k. 

Let  N  =  kv,  where  v  is  the  sample  size  from  each  of  k  populations.  Let  p(u)  = 
(pi(v),p2(v):. . .  ,Pk(v)),  with  Pj(v)  defined  as  in  equation  (2.9).  Then  under  our 
model  with  Fj"'1  — »  F 

— >  avc  f  kS  1 

PCS  =  lim  Pr{pi(u)  >  pj(v)}  =  Pr  <  max  Qt  <  — ^  —  /(2.13) 

^°°  V=2 ••••«*  yj2/ (2k  —  1)  jV 

where 


<?2 

(o\ 

(  1  1/2  ^ 

~  MVN 

1 

\  Qk  ) 

^  1/2  1  ) 

Proof:  See  Appendix  A. 

2.5.4  Combining  BEM  and  AVC  Results 

— >  avc  — ►  bem 

Let  PCS  and  PCS  represent  the  asymptotic  PCS  for  AVC  and  BEM  respec¬ 
tively,  under  the  setup  described  in  §2.5.2  and  §2.5.3.  Combining  the  results  from 
Equations  (2.10)  and  (2.13)  we  have 

- — >  avc  — ►  bem 

PCS  >  PCS 


since 


for  k  >  2  with  6AVC  —  Sbem. 


k8 avc  ^  kShem 

FFF  - 1)  “  sJFk 

Then  we  have 

— y  avc  — y  bem 

PCS  =  PCS 


(2.14) 
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if  and  only  if 


Mavc 


Solving  for  5avc  we  have 


l/2/(2fc  - 1) 


v^T*' 


<Sa“  =  Sbm^k/(2k  -  1). 


(2.15) 


We  use  5bem  and  #bem  (5avc  and  #avc)  to  represent  the  difference,  —  p[fc_1],  or  the 
ratio,  P[jb]/p[A;-i],  respectively  associated  with  BEM  (AVC)  calculations.  Our  goal  is 
to  use  the  relationship  between  5avc  and  5bem  to  define  a  relationship  between  0avc 
and  8hem  that  also  guarantees  (2.14).  This  will  allow  us  to  use  BEM  calculations 
to  approximate  AVC  results. 

Under  the  LFC  for  BEM,  Equation  (2.1),  we  have 


P[k\  ~  P[k- 1] 


8*  1 

8*  +  k  -  1  8*  +  k  -  1 

8*  -  1 

8*  +  k  -  1  ‘ 


(2.16) 


To  relate  this  difference  in  Equation  (2.16)  to  our  asymptotic  analysis,  we  define 
the  sample-size  standardized  difference  as 


ton  Vv(P[k](v)  -P[*-i](v))  =  (k-  1)5  -  (-5) 

=  kS. 


Since  (2.16)  and  (2.17)  represent  the  same  difference,  we  set 


and  solving  for  8* 


kS  = 


8*  -  1 
8*  +  k  —  1 


l  +  (k-l)kS 
1  -k6 


(2.17) 


(2.18) 


(2.19) 
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Equations  (2.18)  and  (2.19)  define  a  relationship  between  a  ratio  and  a  differ¬ 
ence.  We  can  proceed  from  here,  using  Equations  (2.15)  through  (2.19),  to  relate 
#avc  to  #bem  such  that  (2.14)  holds.  Summarizing  the  required  steps  we  arrive  at 
the  following  algorithm. 

Algorithm  2.1  (LSA  for  6avc) 

1.  Express  k6avc  in  terms  of  9avc  using  (2.18). 

2.  Rewrite  k8avc  from  step  1  in  terms  of  k5bem  using  (2.15). 

3.  Express  9bem  in  terms  of  k5bem  from  step  2  using  (2.19). 

Consider  the  following  illustration.  We  have  a  problem  with  a  specified 
(i 9 *,  k,  P*)  where  we  want  to  find  the  minimum  sample  size  required  with  AVC. 
To  approximate  the  required  sample  size,  v,  for  AVC  we  have 

Vavc(9avc  =  9\  k,  P*)  «  vbem(ebem ,  k,  P*)  (2.20) 

where  vavc  and  vbem  denote  the  v  required  for  AVC  or  BEM  respectively. 

We  can  make  this  approximation  since 

PCSbem  ( ybem{9bem,k,P *))  >  P* 

and  from  (2.14)  we  have 

— >  avc  — >  bem 

PCS  (9avc,  k)  =  PCS  ( 9bem,k ), 

which  leads  us  to 

PCSavc  (u6em(06e"\  k,  P*))  «  P*. 
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Continuing  with  our  example  we  step  through  Algorithm  2.1.  At  step  1  we 


kSavc  = 


Qavc  _  j 

Qavc  +  jfc  _  1‘ 


At  step  2  we  have 


9avc  -  1  \  2k  -  1 
6>a?,c  +  k-l)\  k 


Then  at  step  3  we  have  our  LSA  for  9av' 


1  _  (  6avc-l  \  / 2fc  — 1 


(2.21) 


It  will  always  be  the  case  that  9bem  >  6*. 

If  we  want  to  estimate  an  equivalent  9avc  for  9bem  =9*,  we  make  the  following 


approximation, 


vbem^ebem  p*)  w  Vavc(9avc,  k,  P*), 


(2.22) 


using  the  same  reasoning  as  for  (2.20).  Proceeding  as  in  Algorithm  2.1,  we  start 
by  expressing  5bem  in  terms  of  6avc  and  end  up  writing  9avc  in  terms  of  k6avc  to 


obtain 


^avc  _ 


1  +  (*  —  !)( 

i  _  (  9he' 

l  ^bem 


\Q bem+k_l  J  y  2k— 1 
Qbem~l  \  /  k 

9hem+k—\  I  V  2fc-l 


(2.23) 


It  will  always  be  the  case  that  9avc  <  9*.  This  form  of  the  approximation  has  little 
practical  use.  However,  it  does  reflect  another  benefit  of  AVC  in  terms  of  a  smaller 
9 ,  indicating  the  ability  of  AVC  to  discriminate  smaller  differences  between  the  best 
and  the  next  best  system  with  the  same  value  of  v  as  BEM  for  9*.  This  advantage 
becomes  important  in  a  case  where  we  need  to  detect  as  small  a  difference  as 


possible  with  a  fixed  number  of  vector  replications. 
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We  provide  conversions  using  both  of  these  approximations  for  some  common 
values  of  9*  in  Table  2.1  and  present  simulation  results  testing  the  robustness  of 
the  approximations  in  §2.7. 


k 

9* 

Qbem  _  0* 

^avc 

Qavc  _  Q * 

0b  em 

2 

1.2 

1.1604 

1.2506 

2.0 

1.7479 

2.3798 

3 

1.2 

1.1526 

1.2633 

2.0 

1.7205 

2.4297 

4 

1.2 

1.1494 

1.2689 

2.0 

1.7124 

2.4390 

5 

1.2 

1.1476 

1.2720 

2.0 

1.7092 

2.4400 

10 

1.2 

1.1443 

1.2778 

2.0 

1.7061 

2.4326 

Table  2.1:  Equivalent  9  Values  using  LSA 


2.6  Empirical  Results 

In  order  to  allow  easy  comparison  with  available  BEM  results,  we  selected  pop¬ 
ulation  distributions  for  our  simulations  that  allowed  us  to  control  the  value  of 
6*.  These  distributions  are  the  exponential,  continuous  uniform,  and  the  Bernoulli 
presented  in  §2.4.  In  addition,  to  consider  a  less  peaked  continuous  distribution 
without  the  restricted  range  of  the  continuous  uniform,  we  looked  at  a  set  of  gamma 
distributions  with  shape  parameter  a  =  3. 

As  in  our  analytical  results,  we  consider  population  distributions  that  belong 
to  the  same  parametric  family.  We  arbitrarily  select  Hi  as  the  best  population  and 


34 


8 

k 

1.2 

2.0 

2 

1.2000 

2.0000 

3 

1.1589 

1.7808 

4 

1.1368 

1.6632 

5 

1.1227 

1.5885 

Table  2.2:  Exponential:  Values  of  /i  in  LFC  with  A  =  1 

the  remaining  populations  are  identically  distributed.  Let  Xj  represent  a  random 
observation  from  iTj]  j  =  1,  2, . . . ,  k.  We  have 

PrjBest  Population  Wins}  =  Pr{Ara  >  max(A2, . . . ,  AT)}- 

We  then  define 

_  PrjAA  >  max(A~2, . . . ,  Xk)} 

~  (1  -  Pr{W  >  max(A'2, . . . ,  Xk)})  /(k  -  1) ' 

By  setting  8  =  9*,  we  can  then  fix  one  or  more  parameters  for  one  of  the  distribu¬ 
tions  and  solve  for  the  remaining  parameter  to  carry  out  our  simulations  at  a  given 
6*.  Tables  2. 2-2. 5  list  parameters  for  9  =  1.2  and  9  =  2.0  with  k  =  2,  3, . . . ,  5  for 
each  of  our  four  distributions. 

Our  simulation  consists  of  the  following  steps. 

1.  Model  all  systems  using  the  same  distribution  family,  with  system  1  arbi¬ 
trarily  the  best,  and  all  remaining  systems  identically  distributed  such  that 
6  =  9*.  Initialize  PCSbem  and  PCSavc  to  0  and  set  v  =  2. 

2.  Generate  a  set  of  v  random  vector  replications,  where  each  replication  con¬ 
tains  one  observation  for  each  of  the  k  systems. 
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e 

k 

1.2 

2.0 

2 

1.1021 

1.4442 

3 

1.0855 

1.3751 

4 

1.0760 

1.3340 

5 

1.0698 

1.3061 

Table  2.3:  Gamma:  Values  of  j3x  in  LFC  with  ax  =  a0  =  3,  (30  =  1 


e 

k 

1.2 

2.0 

2 

1.1000 

1.5000 

3 

1.0667 

1.3333 

4 

1.0500 

1.2500 

5 

1.0400 

1.2000 

Table  2.4:  Continuous  Uniform:  Values  of  B  in  LFC  with  A  =  1 


e 

k 

1.2 

2.0 

2 

0.4091 

0.1667 

3 

0.4208 

0.2192 

4 

0.4250 

0.2426 

5 

0.4264 

0.2545 

Table  2.5:  Bernoulli:  Values  of  p0  in  LFC  with  px  =  .5 
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3.  For  BEM,  group  the  observations  across  systems  by  vector  replication  and 
count  up  the  number  of  wins  for  each  system.  These  are  our  Yj,  j  = 
1,2,...,*. 

4.  For  AVC,  form  the  vk  pseudo-replications  from  the  v  vector  replications  and 
count  the  number  of  wins  for  each  system.  These  are  our  Zj,  j  =  1, 2, . . . ,  k. 

5.  If  Yi  (BEM  count  associated  with  the  best  system)  is  larger  than  Yj,  j  — 

2, 3 increase  PCSbem  by  1.  If  Y\  ties  for  the  largest  count  with  t  other 
systems,  t  =  1, 2, . . . ,  k  —  1,  increase  PCSbem  by  l/(t  +  1).  If  Y\  <  Yj,  for  any 
j,  j  =  2, 3,  do  not  increase  PCSbem. 

6.  If  Z\  (AVC  count  associated  with  the  best  system)  is  larger  than  Zj,  j  — 
2,3,...,  k,  increase  PCSavc  by  1.  If  Zy  ties  for  the  largest  count  with  t  other 
systems,  t  =  1, 2, . . . ,  k  —  1,  increase  PCSavc  by  l/(t  +  1).  If  Z\  <  Zj,  for  any 
j,  j  =  2,3, ...  ,k,  do  not  increase  PCSavc. 

7.  Repeat  steps  2-6  for  M  macro-replications.  Compute  PCSbem  =  PCSbem/M 
and  PCSavc  =  PCS avc/M. 

8.  Increase  v  and  repeat  steps  2-7. 

Using  parameter  values  from  Tables  2. 2-2. 5,  we  estimated  PCSbem  and  PCSavc 
using  the  simulation  described  above  for  k  =  2, 3, . . . ,  5  populations  out  to  v  =  50 
vectors  for  each  of  the  three  continuous  distributions  at  6  —  1.2  and  2.0.  Results 
were  also  obtained  for  k  =  10  for  exponential  populations.  Due  to  limited  com¬ 
puter  time,  Bernoulli  distributions  were  only  simulated  for  *  =  2  and  3  populations 
at  6  =  1.2.  All  simulation  results  are  for  M  =  100,  000  macro- replications  using  a 
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separate  random  number  stream  for  each  population,  but  common  random  num¬ 
bers  across  distributions.  Since  the  exponential  and  uniform  random  variates  each 
require  a  single  random  number,  results  for  these  distributions  have  synchronized 
random  variates  for  all  of  the  data  generated.  This  is  not  the  case  for  the  gamma 
random  variates  since  an  acceptance-rejection  method  was  used  (Law  and  Kelton 
1991).  The  Bernoulli  random  variates  are  synchronized;  however,  because  of  the 
possibility  of  ties  among  the  systems  within  a  vector  replication,  an  additional 
random  number  stream  is  used  a  random  number  of  times  to  break  the  ties  for  the 
PCS  calculations.  Standard  errors  for  the  PCS  values  are  on  the  order  of  0.0015. 
More  complete  results  are  available  in  Appendix  B. 

Table  2.6  list  results  for  each  of  our  distributions  out  to  v  =  50  vectors  for 
k  =  2  populations  at  9  =  1.2.  The  PCSbem  column  is  from  simulations  using 
exponential  populations.  The  difference  in  the  PCSavc  values  among  the  continuous 
distributions  is  generally  found  in  the  third  decimal  place.  However,  we  see  a 
more  significant  difference  between  the  Bernoulli  PCSavc  and  any  of  the  continuous 
PCSavc  values.  Figure  2.4  demonstrates  the  distributional  dependence  of  PCSavc 
for  exponential  and  Bernoulli  populations.  We  also  notice  significant  improvement 
in  Table  2.6  with  PCSavc  over  PCSbem  for  all  of  the  distributions.  Figure  2.5 
illustrates  the  improvement  with  PCSavc  over  PCSbem  for  k  =  2  to  5  exponential 
populations.  Looking  closely  at  Figure  2.5,  the  spread  between  PCSavc  and  PCSbem 
appears  to  be  increasing  slightly  as  k  increases.  This  is  most  readily  apparent  when 
comparing  the  k  =  2  results  to  the  k  =  3  results.  It  is  also  apparent  from  both 
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PCSavc 

V 

PCSbem 

Exponential 

Uniform 

Gamma 

Bernoulli 

2 

0.5430 

0.5532 

0.5579 

0.5532 

0.5565 

3 

0.5641 

0.5667 

0.5712 

0.5667 

0.5790 

4 

0.5656 

0.5779 

0.5821 

0.5762 

0.5919 

5 

0.5826 

0.5873 

0.5913 

0.5853 

0.6037 

6 

0.5835 

0.5950 

0.5990 

0.5940 

0.6148 

7 

0.5971 

0.6032 

0.6074 

0.6019 

0.6251 

8 

0.5973 

0.6101 

0.6140 

0.6098 

0.6356 

9 

0.6102 

0.6172 

0.6201 

0.6164 

0.6437 

10 

0.6099 

0.6232 

0.6261 

0.6238 

0.6513 

11 

0.6190 

0.6297 

0.6317 

0.6295 

0.6586 

12 

0.6193 

0.6351 

0.6364 

0.6351 

0.6660 

13 

0.6287 

0.6406 

0.6422 

0.6409 

0.6726 

14 

0.6290 

0.6463 

0.6470 

0.6460 

0.6781 

15 

0.6372 

0.6507 

0.6523 

0.6518 

0.6856 

16 

0.6380 

0.6553 

0.6563 

0.6564 

0.6900 

17 

0.6460 

0.6598 

0.6611 

0.6612 

0.6970 

18 

0.6462 

0.6638 

0.6653 

0.6662 

0.7034 

19 

0.6548 

0.6681 

0.6700 

0.6700 

0.7085 

20 

0.6541 

0.6730 

0.6744 

0.6731 

0.7132 

21 

0.6620 

0.6773 

0.6780 

0.6776 

0.7180 

22 

0.6624 

0.6812 

0.6822 

0.6813 

0.7223 

23 

0.6705 

0.6848 

0.6862 

0.6851 

0.7275 

24 

0.6695 

0.6880 

0.6902 

0.6888 

0.7320 

25 

0.6760 

0.6915 

0.6940 

0.6924 

0.7363 

26 

0.6757 

0.6950 

0.6976 

0.6957 

0.7407 

27 

0.6818 

0.6992 

0.7011 

0.6995 

0.7454 

28 

0.6818 

0.7020 

0.7043 

0.7030 

0.7502 

29 

0.6882 

0.7049 

0.7073 

0.7055 

0.7548 

30 

0.6882 

0.7078 

0.7104 

0.7088 

0.7576 

31 

0.6935 

0.7110 

0.7138 

0.7117 

0.7613 

32 

0.6940 

0.7145 

0.7171 

0.7144 

0.7642 

33 

0.6992 

0.7176 

0.7199 

0.7175 

0.7675 

34 

0.6987 

0.7206 

0.7228 

0.7200 

0.7716 

35 

0.7049 

0.7235 

0.7257 

0.7233 

0.7746 

36 

0.7044 

0.7261 

0.7288 

0.7258 

0.7785 

37 

0.7096 

0.7295 

0.7314 

0.7284 

0.7818 

38 

0.7094 

0.7324 

0.7346 

0.7314 

0.7853 

39 

0.7140 

0.7352 

0.7377 

0.7344 

0.7880 

40 

0.7138 

0.7378 

0.7403 

0.7374 

0.7907 

41 

0.7195 

0.7406 

0.7433 

0.7401 

0.7928 

42 

0.7198 

0.7427 

0.7460 

0.7429 

0.7955 

43 

0.7246 

0.7455 

0.7483 

0.7457 

0.7988 

44 

0.7241 

0.7479 

0.7508 

0.7480 

0.8015 

45 

0.7282 

0.7503 

0.7528 

0.7506 

0.8040 

46 

0.7279 

0.7528 

0.7549 

0.7529 

0.8073 

47 

0.7319 

0.7550 

0.7570 

0.7551 

0.8098 

48 

0.7319 

0.7567 

0.7595 

0.7572 

0.8121 

49 

0.7357 

0.7591 

0.7618 

0.7591 

0.8155 

50 

0.7362 

0.7614 

0.7638 

0.7612 

0.8179 

Table  2.6:  PCS  Results  for  k  =  2  Populations  with  6  =  1.2 
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Figure  2.4:  PCSavc  for  Exponential  and  Bernoulli  Populations:  k  =  2  ,9  =  1.2 

Figures  2.4  and  2.5  that  the  spread  between  the  PCSavc  and  PCSbem  values  widens 
as  v  increases.  However,  we  know  that  as  v  approaches  infinity  both  PCSavc  and 
PCSbem  approach  1,  so  this  spread  will  eventually  go  to  zero. 

These  results  clearly  show  an  improvement  in  PCS  with  AVC  for  all  values 
of  k  and  v  considered,  and  also  illustrate  the  weak  dependence  of  PCSavc  on  the 
underlying  population  distributions. 


2.7  Robustness  of  LSA 

To  check  the  accuracy  of  our  LSA  approximation,  we  performed  a  simulation  study. 
The  study  covers  a  range  of  values  for  P*  (.75,  .90,  and  .95)  and  9*  (1.2  and  2.0) 
with  exponential,  continuous  uniform  and  gamma  distributions  for  k  =  2, 3, . . . ,  5 
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Figure  2.5:  PCS  for  Exponential  Populations,  9  =  1.2 

populations  and  Bernoulli  distributions  for  k  =  2,3.  For  the  exponential  and 
continuous  uniform  distributions,  results  are  also  presented  for  k  =  10  at  9*  =  2.0. 
Results  are  not  available  for  all  distributions  at  9*  =  1.2  for  k  >  2.  This  is  due 
to  the  significant  amount  of  computing  time  required  to  obtain  these  results  due 
to  the  much  larger  number  of  replications  required  than  for  9*  =  2.0.  We  have 
included  all  results  available  for  9*  =  1.2  in  Tables  2.7  and  2.9. 

We  first  consider  the  approximation  in  (2.20)  and  perform  the  following  steps. 

1.  Select  a  k  and  9*  and  set  9avc  =  9*.  This  indicates  that  we  are  interested  in 
calculating  AVC  results  at  9*. 

2.  Solve  for  9bem  using  Equation  (2.21). 
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3.  The  calculated  value  of  9bem  will  not  be  in  a  standard  BEM  table.  Use 
FORTRAN  code  developed  by  Goldsman  (1995)  to  find  vbem(9bem,k,P*)  for 
P*  =  .75,  .90,  .95.  Denote  these  values  as  vbem(.75),  vbem(.90),  and  u6em(.95), 
respectively. 

4.  Perforin  simulation  runs  to  estimate  PCSavc  values  at  k  and  9avc  —  9*  when 
using  the  vbem  values  from  step  3.  We  are  looking  for  the  following: 

PCSavc  (vbem{.75)- (9avc,k))  «  .75 
PCSavc  (vbem(.90);  (9avc,  k))  «  .90 
PCSavc  (vbem{.95);  (9avc,k))  «  .95. 

Values  estimated  in  step  4  above  are  reported  in  Tables  2.7  and  2.8.  If  our  LSA 
is  good,  all  PCSavc  values  in  Tables  2.7  and  2.8  should  be  close  to  the  P*  listed  at 
the  top  of  the  column  in  which  they  appear.  The  table  values  include  the  PCSavc 
value  and  the  associated  standard  error  in  parentheses.  All  simulation  runs  use 
the  model  described  in  §2.6  and  are  based  on  M=100,000  macro-replications  for 
the  three  values  of  v  found  in  step  3  above. 

To  illustrate  how  this  approximation  works  for  a  numerical  example,  say 
we  want  to  find  vavc(Qavc  —  1.2,  A;  =  3,  P*  =  .90).  Using  (2.21)  we  obtain 
9bem  =  1.2633,  and  using  FORTRAN  code  developed  by  Goldsman  (1995),  we  find 
vbem(Qbem  _  i .2633,  k  =  3,  P*  =  .90)  =  264.  To  show  how  good  an  approximation 
this  provides  for  our  specified  vavc,  we  simulate  M  =  100,000  macro-replications 
each  containing  264  vector  replications  using  9avc  =  1.2  with  exponential  popu¬ 
lations,  and  obtain  PCSavC  =  .8997  with  a  standard  error  of  0.0009.  So  here  our 
LSA  is  very  good.  These  results  are  included  in  Table  2.7,  where  we  see  that  the 
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k 

^bem 

Distribution 

P*  =  .75 

P*  =  .90 

P*  =  .95 

2 

1.2506 

Exponential 

.7502  (.0014) 

.9012  (.0009) 

.9510  (.0007) 

2 

1.2506 

Gamma 

.7481  (.0014) 

.8994  (.0010) 

.9493  (.0007) 

2 

1.2506 

Uniform 

.7499  (.0014) 

.9010  (.0009) 

.9506  (.0007) 

2 

1.2506 

Bernoulli 

.7856  (.0013) 

.9318  (.0008) 

.9723  (.0005) 

3 

1.2633 

Exponential 

.7470  (.0014) 

.8997  (.0009) 

3 

1.2633 

Gamma 

.7446  (.0014) 

.8994  (.0010) 

■WimM 

3 

1.2633 

Uniform 

.7474  (.0014) 

.8989  (.0010) 

IE 

4 

1.2689 

Exponential 

.7469  (.0014) 

.8964  (.0010) 

.9485  (.0007) 

4 

1.2689 

Gamma 

.7462  (.0014) 

.8984  (.0010) 

.9489  (.0007) 

4 

1.2689 

Uniform 

.7469  (.0014) 

.8948  (.0010) 

.9467  (.0007) 

Table  2.7:  PCS  Achieved  for  6avc  —  9*  =  1.2  Using  LSA 


k 

^bem 

Distribution 

P*  =  .75 

P*  =  .90 

P*  =  .95 

2 

2.3798 

Exponential 

.7614  (.0013) 

.8914  (.0010) 

.9450  (.0007) 

2 

2.3798 

Gamma 

.7608  (.0013) 

.8924  (.0010) 

.9456  (.0007) 

2 

2.3798 

Uniform 

.7625  (.0013) 

.8889  (.0010) 

.9426  (.0007) 

2 

2.3798 

Bernoulli 

.7673  (.0013) 

.9380  (.0008) 

.9773  (.0005) 

3 

2.4297 

Exponential 

.7390  (.0014) 

.8909  (.0010) 

.9426  (.0007) 

3 

2.4297 

Gamma 

.7380  (.0014) 

.8910  (.0010) 

.9440  (.0007) 

3 

2.4297 

Uniform 

.7391  (.0014) 

.8873  (.0010) 

.9391  (.0007) 

3 

2.4297 

Bernoulli 

.7833  (.0013) 

.9279  (.0008) 

.9687  (.0007) 

4 

2.4390 

Exponential 

.7461  (.0014) 

.8911  (.0010) 

.9464  (.0007) 

4 

2.4390 

Gamma 

.7457  (.0014) 

.8929  (.0010) 

.9391  (.0008) 

4 

2.4390 

Uniform 

.7446  (.0014) 

.8858  (.0010) 

.9409  (.0007) 

5 

2.4400 

Exponential 

.7485  (.0014) 

.8910  (.0010) 

.9443  (.0007) 

5 

2.4400 

Gamma 

.7473  (.0014) 

.8940  (.0010) 

.9454  (.0007) 

5 

2.4400 

Uniform 

.7454  (.0014) 

.8824  (.0010) 

.9372  (.0007) 

Table  2.8:  PCS  Achieved  for  6avc  =  9*  =  2.0  Using  LSA 
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results  for  all  the  distributions  achieve  the  desired  P*  to  the  second  decimal  place 
in  almost  all  cases.  In  fact,  we  note  that  the  Bernoulli  results  are  significantly 
larger  than  P*  in  many  cases.  The  cases  where  we  see  more  of  a  departure  from 
P*  are  for  smaller  values  of  P*  where  v  is  typically  less  than  30. 

We  also  notice  that  vbem(Qbem  =  1.2,  k  =  3,  P*  —  .90)  =  437  (Bechhofer,Santner, 
and  Goldsman  1995).  Comparing  this  with  our  approximate  vavc(9avc  =  1.2,  k  = 
3,  P*  =  .90)  =  264,  we  see  a  nearly  40%  reduction  in  the  number  of  replications 
required  with  AVC. 

To  be  complete  we  consider  using  the  approximation  in  (2.22)  and  perform  the 
following  steps. 

1.  Select  a  k  and  9*  and  set  9bem  =  9*.  This  indicates  that  we  are  interested  in 
using  BEM  results  at  9*. 

2.  Solve  for  9avc  using  Equation  (2.23). 

3.  Go  to  standard  BEM  tables,  such  as  Table  8.1  in  Bechhofer,  Santner, 
and  Goldsman  (1995),  and  read  off  values  for  vbem (9bem ,  k ,  P*)  for  P*  = 
.75,  .90,  .95.  Denote  these  values  as  vbem(.75),  vbem(.90),  and  vbem(.95)  re¬ 
spectively. 

4.  Perform  simulation  runs  to  obtain  PCSavc  values  at  k  and  9avc  from  step  2 
when  using  the  vbem  values  from  step  3.  We  are  looking  for  the  following: 

PCSavc  (vbem(.75);  ( 9avc ,  k))  «  .75 
PCSavc  (vbem(.90);  ( 9avc ,  k))  ^  .90 
PCSavc  (u6em(.95);  ( 9avc ,  k))  «  .95. 
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Values  estimated  in  step  4  above  are  reported  in  Tables  2.9  and  2.10.  The  table 
values  include  the  PCSavc  value  and  the  associated  standard  error  in  parentheses. 
All  simulation  runs  use  the  model  described  in  §2.6  and  are  based  on  M=100,000 
macro-replications  for  the  three  values  of  v  found  in  step  3  above. 

Following  our  previous  numerical  illustration,  suppose  we  want  to  find 
vbem^0bem  _  ^  2,  k  =  3,  P*  =  .90).  Using  (2.23)  we  obtain  9avc  =  1.1526.  We 
go  to  a  standard  BEM  table,  to  find  our  specified  vbem  =  437.  As  we  did  above,  we 
test  the  robustness  of  our  approximation  by  simulating  100,000  macro-replications 
each  containing  437  vector  replications  using  9avc  =  1.1526  with  exponential  pop¬ 
ulations,  and  obtain  PCSavc  =  .8991  with  a  standard  error  of  0.0010.  Once  again 
our  LSA  is  very  good.  These  results  are  presented  with  similar  results  for  our  other 
selected  distributions  in  Tables  2.9  and  2.10. 

The  benefit  from  this  form  of  the  LSA  is  reflected  by  9avc  <  9*.  This  indicates 
that  AVC  can  provide  better  discrimination  between  the  systems  at  the  same  level 
of  confidence  and  with  the  same  data. 

2.8  Conclusions 

When  trying  to  pick  the  best  system  out  of  k  systems,  there  are  many  instances 
when  this  selection  should  be  based  on  one-time  performance  rather  than  long-run 
average  performance.  Multinomial  selection  procedures  provide  a  framework  for 
defining  such  a  problem,  and  Procedure  BEM  is  the  classical  approach  for  solving 
it.  Procedure  AVC  is  an  alternative  approach  designed  to  obtain  a  higher  PCS 
by  performing  all  possible  comparisons  across  all  systems  for  a  given  set  of  system 
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k 

Q&vc 

Distribution 

P*  =  .75 

P*  =  .90 

P*  =  .95 

2 

1.1604 

Exponential 

.7486  (.0014) 

.8998  (.0009) 

.9500  (.0007) 

2 

1.1604 

Gamma 

.7493  (.0014) 

.9008  (.0009) 

.9495  (.0007) 

2 

1.1604 

Uniform 

.7505  (.0013) 

.9011  (.0009) 

.9500  (.0007) 

3 

1.1526 

Exponential 

.7514  (.0014) 

.8991  (.0010) 

.9500  (.0007) 

3 

1.1526 

Gamma 

.7502  (.0014) 

.8997  (.0010) 

.9500  (.0007) 

3 

1.1526 

Uniform 

.7510  (.0013) 

.8983  (.0009) 

.9493  (.0007) 

4 

1.1494 

Exponential 

.7504  (.0014) 

.8998  (.0010) 

.9501  (.0007) 

4 

1.1494 

Gamma 

.7520  (.0014) 

.9002  (.0009) 

.9509  (.0007) 

4 

1.1494 

Uniform 

.7499  (.0014) 

.8975  (.0010) 

.9489  (.0007) 

Table  2.9:  PCS  Achieved  for  6bem  =  9*  =  1.2  Using  LSA 


k 

^avc 

Distribution 

P*  =  .75 

P*  =  .90 

P*  =  .95  • 

2 

1.7479 

Exponential 

.7692  (.0013) 

.9021  (.0009) 

.9471  (.0007) 

2 

1.7479 

Gamma 

.7702  (.0013) 

.9037  (.0009) 

.9477  (.0007) 

2 

1.7479 

Uniform 

.7703  (.0013) 

.9013  (.0009) 

.9453  (.0007) 

2 

1.7479 

Bernoulli 

.7950  (.0013) 

.9373  (.0008) 

.9726  (.0005) 

3 

1.7205 

Exponential 

.7471  (.0014) 

.8971  (.0010) 

.9460  (.0007) 

3 

1.7205 

Gamma 

.7456  (.0013) 

.8975  (.0009) 

.9464  (.0007) 

3 

1.7205 

Uniform 

.7459  (.0014) 

.8933  (.0010) 

.9432  (.0007) 

3 

1.7205 

Bernoulli 

.7911  (.0013) 

.9332  (.0008) 

.9701  (.0005) 

4 

1.7124 

Exponential 

.7419  (.0014) 

.8944  (.0010) 

.9461  (.0007) 

4 

1.7124 

Gamma 

.7418  (.0014) 

.8956  (.0010) 

.9464  (.0007) 

4 

1.7124 

Uniform 

.7403  (.0014) 

.8888  (.0010) 

.9414  (.0007) 

5 

1.7092 

Exponential 

.7458  (.0014) 

.8950  (.0010) 

.9474  (.0007) 

5 

1.7092 

Gamma 

.7452  (.0014) 

.8976  (.0010) 

.9478  (.0007) 

5 

1.7092 

Uniform 

.7412  (.0014) 

.8881  (.0010) 

.9410  (.0007) 

10 

1.7061 

Exponential 

.7576  (.0014) 

.9024  (.0009) 

.9516  (.0007) 

10 

1.7061 

Uniform 

.7476  (.0014) 

.8891  (.0009) 

.9411  (.0007) 

Table  2.10:  PCS  Achieved  for  0bem  —  9*  =  2.0  Using  LSA 
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performance  data.  Construction  of  Procedure  AVC  closely  follows  that  of  BEM, 
allowing  researchers  to  easily  move  from  a  standard  approach  to  our  new  approach. 


p* 

k 

9* 

.75 

.90 

.95 

1.01 

12171  (18371) 

* 

* 

1.05 

509  (765) 

1839  (2759) 

3027  (4545) 

2 

1.10 

133  (201) 

483  (723) 

793  (1191) 

1.20 

37  (55) 

133  (199) 

217  (327) 

2.00 

3(5) 

9(15) 

15  (23) 

1.05 

1544  (2565) 

3741(6211) 

5526(9165) 

3 

1.10 

401  (666) 

972  (1615) 

1436  (2385) 

1.20 

108  (181) 

264  (437) 

388  (645) 

2.00 

7(12) 

17  (29) 

25  (42) 

4 

1.20 

187  (326) 

398  (692) 

565  (979) 

2.00 

12  (20) 

25  (43) 

36  (61) 

5 

vm 

271  (486) 

541  (964) 

748  (1331) 

E9 

17  (29) 

33  (58) 

46  (81) 

Table  2.11:  Minimum  Number  of  Vectors  to  Achieve  P*  for  AVC  (BEM) 


From  the  simulation  design  point  of  view,  AVC  can  also  be  used  to  our  advan¬ 
tage  by  allowing  us  to  plan  a  smaller  number  of  replications  to  achieve  a  desired 
PCS,  P*.  Table  2.8  presents  comparisons  of  the  minimum  number  of  independent 
replications  needed  to  achieve  a  given  P*  for  AVC  and  BEM.  Values  for  BEM 
are  taken  from  Bechhofer,  Santner,  and  Goldsman  (1995).  The  AVC  values  are 
obtained  using  our  LSA  in  (2.21)  with  9avc  =  9*  to  find  9bem  and  then  running 
an  exact  code  for  PCSbem  provided  by  Goldsman  (1995)  at  9  =  9bem.  An  asterik 
indicates  that  runs  were  not  accomplished  due  to  the  large  number  of  replications 
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(>  40,  000)  required.  As  k  increases,  we  see  a  more  dramatic  reduction  in  the 
number  of  vector  observations  needed  with  AVC  to  achieve  the  same  P*.  So  the 
advantages  of  AVC  over  BEM  appear  greater  for  more  challenging  MSPs. 
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CHAPTER  3 


POINT  ESTIMATION  FOR  MULTINOMIAL 
SUCCESS  PROBABILITIES 


3.1  Introduction 

Suppose  we  have  k  >  2  independent  populations,  denoted  Tiy,  7r2, . . . ,  7 t*.  In  a 
simulation  context  each  population  is  a  simulated  system.  We  consider  the  problem 
of  estimating  the  probabilities  associated  with  each  system  being  the  best  system 
in  a  single  comparison  (trial)  of  simulated  results  for  all  of  the  systems. 

Let  Xji  represent  the  ith  replication  from  system  j  of  some  performance  mea¬ 
sure.  Each  system  (tcj,  j  =  1, . . . ,  k)  has  an  unknown  constant  probability  (pj,  j  = 
l, ...  ,k)  of  having  the  largest  value  of  the  performance  measure  in  any  replication 
containing  one  observation  from  each  system.  We  define  the  best  system  as  the 
system  most  likely  to  have  the  largest  performance  measure  (i.e.,  it  wins)  in  any 
comparison  across  all  systems.  Such  a  comparison  corresponds  to  a  multinomial 
trial,  where  one  and  only  one  system  can  win  in  any  given  trial  and  the  numbers 
of  wins  for  the  individual  systems  in  v  independent  trials  follows  a  multinomial 
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distribution.  Our  objective  is  to  provide  estimates  for  the  unknown  multinomial 
success  probabilites,  pj,  j  =  1, . . . ,  k. 

Applications  where  a  multinomial  distribution  is  appropriate  include  selecting 
the  best  of  a  set  of  tactical  or  strategic  military  actions,  where  best  means  the  max¬ 
imum  damage  in  a  single  strike.  An  example  in  the  area  of  structural  engineering 
is  finding  the  design  that  performs  best  in  a  one-time  catastrophic  event,  such  as 
an  earthquake.  Simulation  examples  include  selecting  the  schedule  most  likely  to 
result  in  completing  all  jobs  on  time;  selecting  the  investment  portfolio  most  likely 
to  provide  the  largest  return;  or  selecting  the  computer  system  with  the  highest 
probability  of  completing  a  series  of  tasks  without  failure.  Each  of  these  appli¬ 
cations  involves  the  comparison  of  quantitative  measures  of  performance  among 
competing  systems  as  opposed  to  comparing  qualitative  measures.  For  the  type  of 
application  considered  in  this  study,  we  require  a  quantitative  measure  of  system 
performance  for  each  system  on  each  trial  to  be  compared  with  the  performance 
of  other  systems  across  any  or  all  of  the  remaining  trials. 

Let  X;  =  (A'  12,  A-22,  • . . ,  A represent  the  1  replication  across  all  k  systems. 
Let  Yji  =  1  if  Xji  >  X^i,  for  i  —  1,  2, . . . ,  k,  but  £  ^  j;  and  let  Yji  =  0  otherwise.  In 
other  words,  Yji  =  1  if  Xji  is  the  largest  observation  in  X.(.  In  case  of  a  tie  for  the 
largest  value,  we  randomly  select  one  of  the  tied  populations  as  the  best.  Suppose 
that  there  are  v  independent  replications  across  all  systems,  and  let  Yj  =  £V=1  Yji 
represent  the  number  of  times  system  j  wins  out  of  these  v  replications.  Let 
Pj  =  Pr{Xjj  >  Xu,  W  ^  ]  }  where  0  <  p;/  <  1  and  Y.kj=\Pj  =  1-  Then  £j=i  Yj  =  v 


51 


and  the  k- variate  discrete  random  variable  Y  =  (T|,  Y2, . . . ,  Yk)  follows  a  multino¬ 
mial  distribution  with  v  trials  and  success  probabilities  p  =  (pi,p2,  ■  ■  ■  ,Pk)-  The 
probability  mass  function  for  Y  with  parameters  v  and  p  is 

Pi'O'i  =  !Ji ,  Y2  =  y2,...,  Yk  =  yk)  =  — ^ — -  J]  pf. 

Hj=i  y.r  j= 1 

One  set  of  point  estimators  for  the  multinomial  success  probabilities  are  the  max¬ 
imum  likelihood  estimates  (MLEs)  given  by  pj  =  Yj/v,  j  =  1,  2, . . . ,  k. 

Closely  related  to  the  estimation  of  the  pjs  is  the  problem  of  determining  which 
of  the  systems  has  the  largest  probability  of  being  the  best  system.  This  is  known 
as  the  multinomial  selection  problem  (MSP).  The  classical  solution  procedure  for 
the  MSP,  Procedure  BEM  (Bechhofer,  Elmaghraby,  and  Morse  1959),  uses  the 
multinomial  cell  counts,  Y,  to  select  the  best  system.  Some  of  the  concepts  and 
terminology  used  in  the  MSP  carry  over  to  our  point  estimation  problem. 

Due  to  convention  and  convenience  when  comparing  simulated  system  re¬ 
sponses,  the  responses  are  typically  grouped  by  replication,  corresponding  to  a 
trial  in  a  physical  experiment.  Grouping  system  responses  in  this  fashion  is  ar¬ 
bitrary  and  since  our  simulated  responses  are  quantitative,  we  can  compare  any 
observation  from  one  system  with  any  observation  from  each  of  the  remaining  sys¬ 
tems.  This  means  that  a  single  observation  from  system  1  can  be  grouped  in  a 
vector  comparison  with  any  one  of  the  v  observations  from  system  2,  and  with  any 
one  of  the  v  observations  from  system  3,  and  so  on,  up  to  and  including  any  one 
of  the  v  observations  from  system  k.  Since  there  are  v  observations  from  system  1 
as  well,  this  gives  us  a  total  of  vk  vector  comparisons  (trials)  that  can  be  formed 
with  v  independent  observations  from  each  of  the  k  systems.  We  incorporate  this 


52 


setup,  which  we  call  AVC  for  All  Vector  Comparisons,  to  construct  new  point  es¬ 
timators  for  the  pjS.  Our  estimators  turn  out  to  be  fc-sample  U-statistics  (Randles 
and  Wolfe  1979).  By  arbitrarily  performing  only  the  v  vector  comparisons  where 
the  observations  for  each  system  are  from  the  same  replication,  we  disregard  the 
information  available  from  the  remaining  vk  —  v  comparisons. 

Our  results  suggest  some  advantages  of  AVC  estimators  over  the  standard 
MLEs.  First,  we  prove  that  the  variance  of  the  AVC  estimators  is  no  larger  than 
the  variance  of  the  MLEs.  For  specific  small-sample  examples  we  demonstrate  the 
magnitude  of  the  variance  reduction  with  AVC  estimators  for  a  fixed  v.  We  show 
this  analytically  for  small  values  of  v  and  We  also  show  that  the  variance  re¬ 
duction  depends  weakly  on  the  distributions  of  the  system  performance  measures. 
This  fact,  along  with  the  difficulty  of  analytically  evaluating  the  exact  variance  for 
AVC  estimators  for  even  small  A;  and  v,  led  us  to  a  consider  using  the  asymptotic 
variance  as  an  approximation  to  the  exact  variance.  An  estimate  of  the  asymptotic 
variance  is  used  to  form  confidence  intervals. 

This  paper  is  organized  as  follows:  We  first  briefly  address  the  MSP  and  intro¬ 
duce  terminology  used  in  our  point  estimation  problem.  We  then  present  the  MLEs 
and  the  associated  variances  for  the  PjS.  Then  we  describe  our  new  approach,  AVC, 
which  represents  our  point  estimators  as  A;-sample  U-statistics.  We  present  a  proof 
for  a  variance  reduction  with  AVC  along  with  analytical  results  showing  the  mag¬ 
nitude  of  the  variance  reduction  for  small  k  and  v  covering  different  population 
distributions  for  the  performance  measures.  Our  asymptotic  approximation  for  the 
variances  of  the  AVC  estimators  is  then  presented.  Empirical  results  follow  and 
include  the  construction  of  confidence  intervals. 
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3.2  Background 


Because  of  the  close  relationship  between  the  MSP  and  the  estimation  of  multino¬ 
mial  success  probabilities,  we  briefly  address  a  classical  solution  procedure  for  the 
MSP  and  define  some  terminology  we  will  be  using  throughout  our  discussion. 

Bechhofer,  Elmaghraby  and  Morse  (1959)  describe  a  single-stage  procedure, 
Procedure  BEM,  for  selecting  the  multinomial  event  (population  or  system)  which 
has  the  largest  success  probability.  Procedure  BEM  requires  the  specification  of 
P*  (where  1  jk  <  P*  <  1),  a  minimum  probability  of  correctly  identifying  the 
population  with  the  largest  success  probability  (i.e.,  the  best  population),  and  9* 
(where  1  <  9*  <  oo),  the  minimum  ratio  of  the  largest  success  probability  to 
the  second  largest  success  probability  that  we  want  to  be  able  to  detect.  The 
probability  of  correct  selection  (PCS)  is  a  property  of  the  selection  procedure  and 
provides  no  information  about  the  values  of  the  pjS.  The  procedure,  as  adapted  to 
simulation,  consists  of  the  following  steps: 

Procedure  3.1  (BEM) 

1.  For  given  k  and  9* ,  find  the  minimum,  value  of  v,  denoted  as  v* ,  that  guar¬ 
antees  that  the  PCS  is  at  least  P*. 

2.  Generate  v*  independent  replications  for  each  population. 

3.  Compute  Yj  =  £i=i  Yji}  for  j  =  1, 2, . . . ,  k. 


54 


4-  Let  l'(i)  <  1(2)  <  ■  •  •  <  Y(k)  be  the  ranked  sample  counts  from  step  3.  Select 
the  population  associated  with  the  largest  count,  Y^k),  as  the  best  population. 
In  case  of  a  tie  for  the  largest  count,  randomly  select  one  of  the  tied  popula¬ 
tions  as  the  best. 

To  determine  the  appropriate  v*  in  step  1,  let  ppj  <  p[2]  <  •  •  •  <  p\k\  denote  the 
ranked  success  probabilities  for  the  k  populations.  Since  only  values  of  the  ratio 
9  —  P[k]/P[k-i)  greater  than  or  equal  to  9*  are  of  interest,  we  are  indifferent  between 
the  best  and  the  next-best  population  for  values  of  6  <  9*.  A  procedure  of  this 
type  is  referred  to  as  an  indifference- zone  approach.  Select  v*  as  the  minimum 
number  of  independent  vector  observations  required  to  achieve  a  PCS  greater  than 
or  equal  to  P*  whenever  9  >9* . 

We  define  the  least  favorable  configuration  (LFC)  of  p  =  (p[i),p[2],  ■  ■  -  ,P[k])  as 
the  configuration  where  PCS  is  a  minimum  over  all  configurations  with  9  >  9* 
(Gibbons,  Olkin,  and  Sobel  1977).  If  we  obtain  a  PCS  >  P*  with  our  selected  v* 
under  the  LFC,  then  a  PCS  of  at  least  P*  can  be  guaranteed  for  any  configuration 
of  p  with  9  >9* .  Keston  and  Morse  (1959)  prove  that  the  LFC  for  BEM  is  given 
by 

1 

pi  i]  =  m  =  ■■■  =  P[k~  i]  =  g,  +  k  _  1 

9* 

P{k]  ~  9*  +  k  -  1 ' 

A  well  known  set  of  estimators  for  our  multinomial  success  probabilities  are  the 
MLEs  given  by 

Pj  =  (3.1) 
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which  represent  the  fraction  of  wins  out  of  v  replications  for  population  j,  j  = 
1 , ...  ,k.  It  is  well  known  that  E[p,]  =  pj  and 

Var(ft)  -  .  (3.2) 

3.3  AVC  Estimators 

We  propose  a  method  to  provide  point  estimators  for  the  multinomial  success 
probabilities  with  smaller  variances  than  the  corresponding  MLEs  using  the  same 
replications  Xj,  i  =  1,2, ...  ,v.  Rather  than  comparing  the  ith  replication  for  each 
system  with  the  ith  replications  of  the  other  systems,  consider  a  total  of  vk  pseudo¬ 
replications  formed  by  associating  each  Xji  (j  =  1, 2, . . . ,  k\  i  =  1, 2, . . . ,  v),  with 
all  possible  combinations  of  the  remaining  AW  (£  =  1,2, ...  ,k]  £  ^  j;  h  = 
1,  2, . . . ,  v).  Each  such  pseudo-replication  contains  one  observation  from  each  pop¬ 
ulation.  Notice  that  the  vk  pseudo-replications  include  the  v  independent  replica¬ 
tions  from  which  the  pseudo- replications  are  formed. 

Define 

n  -  X:,)  (3-3) 

ai  =  la2=l  a*  =  1 

for  j  =  1,  2, . . . ,  k  with 

'1,  a  >  0 

x  0,  a  <  0 

na)  = 

v  randomly  assign 

,  0  or  1,  a  =  0. 

Thus,  Zj  represents  the  number  of  times  out  of  vk  pseudo-replications  that  popu¬ 
lation  7 Tj  wins  (ties  broken  randomly)  and  Ej=i  Zj  —  vk. 
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As  a  specific  illustration  of  how  these  pseudo-replications  are  formed,  consider 
k  =  3  systems  with  two  observations  generated  for  each.  Our  original  replications 


and  counts  (Y),  j  —  1,2,3)  are  then 

An  A21  A31 

Al2  A22  A  32 

IT  A>  aT 

We  will  have  a  total  of  vk  (23  =  8)  pseudo-replications  including  our  two  original 
replications  when  using  AVC.  These  pseudo-replications  and  the  associated  counts 
(Zj,  j  =  1,  2,  3)  for  each  system  are 


An 

a21 

Ar3i 

An 

a21 

A32 

An 

A22 

A3i 

An 

A22 

A32 

A 12 

A21 

A31 

A12 

A21 

A32 

A12 

A22 

A31 

Ai2 

A22 

A32 

Z\  Z‘i  A3 

Our  new  point  estimators  are  then 

=  >  (3-4) 

which  represent  the  fraction  of  wins  out  of  vk  pseudo-replications  for  population 
j.  We  refer  to  these  estimators  in  (3.4)  as  AVC  estimators,  or  to  p  =  (pi, . . .  ,pk) 
collectively  as  our  AVC  estimator.  Clearly  E [pj]  —  pj,  but  the  Var (pj)  is  more 
complex  to  calculate  than  Var (p,).  To  find  the  variances  for  the  individual's,  we 
represent  our  AVC  estimator  as  a  A;-sample  U-statistic,  where  p  =  (pi,P2,  ■  ■  ■  ,Pk ) 
represents  our  parameters  for  the  system  success  probabilities  and  X )  ~  Fj .  From 
Randles  and  Wolfe  (1979,  p.  104),  we  say  p  is  estimable  of  degree  (!,...,!)  for 
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distributions  (Fi, . . . ,  Fk)  in  some  family  F,  if  (1, . . . ,  1)  are  the  smallest  sample 
sizes  from  across  all  systems  or  populations  for  which  there  exists  an  unbiased 
estimator  of  p  for  every  (iq, . . . ,  Fk)  E  F.  Formally  stated  we  have 

for  j  =  1, . . . ,  k  where  we  define  our  &:-sample  symmetric  kernel  h^\-)  as 

k 

hu)  =  n  4>{Xji  —  Xu).  (3.5) 

<=i;  l  ft 

So  we  have  —  1  if  the  observation  from  the  jth  system  is  the  largest  in  any  vector 
comparison  across  all  systems.  Substituting  our  kernel  from  (3.5)  into  (3.3)  and 
dividing  by  vk,  we  have  presented  our  pj,  j  =  1, 2, . . . ,  k  as  /c-sample  U-statistics. 

Using  well-known  results  for  U-statistics  (Randles  and  Wolfe  1979),  we  develop 
the  variance  of  a  2-sample  U-statistic  in  general.  First  we  have 

E(Fi,f2)  [h(Xu,  X2i)]  =  p. 

Then  the  covariance  terms  are 

6,o  =  Cov  [h(Xicn,  X2q2),  h(XlQl)  X2p2)} 

=  E  [/l(AlQl ,  X2a2)h{Xlai,  Xlfo)}  —  p2 

6.1  =  Cov  [h(Xiai ,  A2a2),  A^2a2)] 

=  E[h(Xiai,X2a2)h(Xip1,X2a.2)\  —  p2 

6.1  =  Cov  [h(Xiai ,  X  2a2)i  h(X lQi,A2a2)] 

E  [h(Xiai ,  X2a2)h(Xiai ,  A2a2)]  P 

where  cq  /  Pi  and  6,0  =  0.  In  our  specific  context,  suppose  we  want  to  estimate 
Pi  =  Pr{AT  >  AY}.  We  use  our  kernel  from  equation  (3.5)  for  Yx)(.)  and  define 

6,0  =  Pr{Abai  >  AY*2;  ATai  >  X2p2}  —  p\  (3.6) 


58 


£o,i  = 

Pr{Xlai  >  X2 Q2;  Xifa  >  X2a.2 }  —  p\ 

(3,7) 

6,i  = 

Pr{AdQl  >  V2Q2 ;  Ado,  >  ALn2 }  —  p\ 

(3.8) 

where  oq,  a2,  A,  fa  €  {1, . 

Then  using  a  general  equation  for  the  variance  of  a  2-sample  U-statistic,  we 
can  express  the  variance  of  our  estimators  as 


3.4  Analytical  Results 

The  following  analytical  study  illustrates  a  number  of  important  properties  of  our 
AVC  estimator.  First,  we  provide  a  general  proof  that  the  variance  of  our  AVC 
estimator  is  no  larger  than  the  variance  of  the  MLE.  Then  we  quantify  this  variance 
reduction  with  AVC  for  specific  cases.  We  also  show  that  the  variance  reduction  of 
AVC  relative  to  MLE  depends  weakly  on  the  underlying  distribution  for  the  Xji. 
In  addition  we  demonstrate  the  difficulty  in  obtaining  analytical  results  for  even 
a  small  number  of  populations  and  observations,  and  thus  provide  motivation  for 
our  asymptotic  variance  approximation. 

3.4.1  Smaller  Variance  of  the  AVC  Estimator 

To  show  that  AVC  provides  a  point  estimator  with  a  smaller  variance  than  MLE, 
suppose  we  have  simulation  results  for  k  systems  with  v  observations  from  each 
system  as  shown. 
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An 

X2i  ■ 

••  Xkl 

AT  2 

X‘22  ‘ 

••  xk2 

ATV 

x2v  ■ 

X  kv 

Fi 

f2  ■ 

••  Fk- 

Let  F\,  F2, . . . ,  Fk  represent  the  empirical  cumulative  distribution  function  (cdf) 
for  each  population  and  suppose  we  wish  to  estimate  p\.  Our  MLE  point  estimator 
is  then 


V\ 


with  variance 

Var(pi)  =  ?l(1 

V 

To  arrive  at  our  AVC  estimator  let 


Pi(F)  =  Pr{Xu.  >  A>,  Vj  +  1,  f 


1, . . . ,  vk\ Fi,  F2,...,  Fk} 


where  i*  represents  the  number  of  the  pseudo- replication  and  Pi(F)  is  the  proba¬ 
bility  system  i  is  the  best  when  the  data  are  distributed  as  Fi:  F2, . . . ,  Fk.  Under 
Fi,  F2, . . . ,  Fk,  each  pseudo-replication  occurs  with  an  equal  probability  of  l/vk. 
We  can  then  write 


Pi(F)  = 


(the  number  of  pseudo-replications  where  population  1  wins) 

vk 


implying  that  the  numerator  is  Z\  from  (3.3).  Therefore,  we  can  define  our  AVC 
estimator  as 

Pi  =Pi(F). 
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Using  standard  definitions  for  conditional  expectation  and  variance  (Casella 
and  Berger  1990,  p.  158)  we  know 

Varfo)  =  Vav[E(p1\F1,F2,...,Fk)]+E[Yav(pl\FuF2,...,Fk^ 

=  Var^F)]  +  E[Var(pj \F\,  F2, Fk)] 

=  VM{pl}  +  E[VaT(p1\Fl,F2,...,Fk)}. 

Since  E[Var(pi|Fi,  F2, . . . ,  F*)]  >  0,  we  have  shown  that 

Var(p!)  <  Varfpi ). 

In  the  following  section  we  quantify  the  size  of  this  reduction  in  the  AVC  estimator 
variance  for  specific  cases. 

3.4.2  Small-Sample  Results 

We  restrict  our  attention  to  continuous  distributions  for  the  A^s,  which  eliminates 
the  possibility  of  ties  among  the  observations.  We  let  n[k]  be  the  best  popula¬ 
tion  and  assume  all  the  remaining  populations,  7T[i], . . . ,  K[k-i\  are  identically  dis¬ 
tributed.  This  setup  gives  us  the  LFC  for  BEM.  We  also  consider  all  population 
distributions  to  belong  to  the  same  parametric  family.  For  illustrative  purposes, 
let  X  represent  an  observation  from  7 r[fc]  and  let  O  represent  an  observation  from 
any  of  the  remaining  inferior  populations. 

First,  consider  X  ~  exp(A)  and  O  ~  exp (p)  and  let  A  <  fx,  where  A  >  0  and 
H>  0  are  exponential  rates.  This  particular  example  was  suggested  by  Goldsman 
(1995)  in  some  early  discussions  concerning  this  research.  For  k  =  2  and  v  =  2, 
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we  have  p[2]  =  Pr{A'  >  0}  =  p/(X  +  p)  and  p[ p  =  Pr{AA  <  0}  —  A/(A  +  p).  We 
arbitrarily  let  k  be  the  best  population  so  pk  —  P[k\-  From  (3.2)  we  have 

Xp, 

2(A  +  m)2' 

To  find  Var(p2),  we  must  first  find  the  covariance  terms.  From  (3.6)  we  have 


£i,o  —  Pr{Ab  >  Op  Xi  >  02}  —  pi 

=  Pr{Ah  >  M}-pl 

POO 

=  /  2/ie-(A+M)m(l  -  e“^m)  dm  -  p\ 

Jo 

2p2  p2 

(A  +  2p)(\  +  p)  (A  +  p)2 
_  A  p2 

(A  +  2p)(X  +  p)2 

where  M  =  maximum(01, 02).  Then  from  (3.7) 
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where  N  =  minimum  (AA ,  Ar2) .  Finally,  from  (3.8) 


6,1  =  Pr{X  l>01;X1>01}-pl 

k _ F2 

A  +  n  (A  +  fj,)2 

(A  +  //)2 ' 

Inserting  these  results  into  (3.9)  we  have 

Var(p2)  =  -  [6,o  +  6,i  +  6,o  +  6,i] 

3  Ap(A2  +  3A//  +  p2) 

4  (2A  +  /i)(A  +  2//) (A  +  p)2 


To  illustrate  the  variance  reduction  achieved  by  the  AVC  estimator,  we  display 
the  ratio  of  Var(p2)  to  Var(p2)  plotted  against  6  =  p2/pi  =  yu/A  in  Figure  3.1. 
These  results  are  for  exponential  populations  with  k  =  2  and  v  =  2.  A  ratio 
less  than  1  indicates  a  variance  reduction  with  AVC.  The  reduction  in  variance 
is  on  the  order  of  roughly  20%  over  a  range  of  6  between  1  and  4.  This  covers 
most  of  the  practical  range  of  0  included  in  standard  tables  for  BEM  and  used 
by  experimentors.  Clearly,  as  6  increases,  both  estimators  approach  1  and  the 
associated  variances  approach  0. 

In  order  to  demonstrate  the  weak  dependence  of  the  variance  reduction  on  the 
distribution  of  the  Xji,  consider  the  following  example  for  k  =  2  and  v  =  2  where 
X  and  O  are  from  continuous  uniform  populations.  Specifically  X  ~  U(0,B)  and 
O  ~  [7(0,  A),  where  0  <  A  <  B.  Then  we  have  p[2]  =  Pr{X  >  0}  =  (2 B  -  A)J2B 
and  pp]  =  Pr{Ar  <  0}  =  A/2B.  From  (3.2)  we  have 


Var(p2) 


1  A{2B  -  A) 
8  B~2 
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e 

Figure  3.1:  Ratio  of  Var(p2)/Var(p2)  for  Exponential  Populations  with  k  =  2, 


We  can  then  find  the  corresponding  variance  of  the  AVC  estimator  in  a  similar 
fashion  as  we  did  for  the  exponential  populations.  From  (3.9)  we  obtain 


Var(p2) 


5A(2R  -  A) 

48R2 


In  Figure  3.2  we  show  the  reduction  in  the  variance  of  the  AVC  estimator  for 
our  uniform  populations  again  as  a  ratio  of  Var(p2)/Var(p2).  The  MLE  estimator 
variance  is  identical  for  both  exponential  and  uniform  distributions.  We  see  a 
flatter  curve  in  Figure  3.2  than  in  Figure  3.1.  In  Figure  3.3  we  illustrate  the 
difference  between  the  exponential  and  uniform  populations  by  plotting  the  ratio 
of  Var(p2)  for  the  uniform  populations  to  Var(p2)  for  the  exponential  populations. 
It  is  interesting  to  notice  that  this  difference  is  neglible  over  most  of  the  practical 
range  of  9. 
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Figure  3.2:  Ratio  of  Var(p2)/Var(p2)  for  Uniform  Populations  with  k  =  2,  v  =  2 

Once  we  have  our  covariance  terms  for  a  given  k ,  we  can  find  the  exact  variance 
for  pj  for  any  v.  For  example  with  k  =  2  using  (3.9)  we  have 

Var(pj(n))  =.  —  [{v  —  1)  (£0,i  +  Ci,o)  +  Ci,i]  •  (3-11) 

As  we  increase  k  we  add  additional  covariance  terms  to  our  variance  calculations. 
For  example  at  k  =  3,  we  can  modify  (3.11)  to  get 

Var {Pj(v))  —  ^3  [(w  —  l)2  (Ci, o,o  +  Co, i,o  +  Co,o,i)  + 

( v  ~  1)  (Ci, i,o  +  Ci, o,o  +  Co,i,i)  +  Cl,l,l] 

where  £c,d;e  are  the  analagous  covariance  terms  for  k  —  3  populations,  with  c,d,eE 
{0, 1}.  These  additional  covariance  terms  become  increasingly  difficult  to  calculate 
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Figure  3.3:  Ratio  of  Uniform  Var(p2)  over  Exponential  Var(p2)  for  k  =  2,  v  =  2 

for  even  moderately  large  k.  This  computational  complexity,  along  with  the  weak 
distributional  dependence  of  the  AVC  variance  reduction,  led  us  to  consider  an 
asymptotic  approximation  for  the  £>sample  U-statistic  variance. 

3.5  Asymptotic  Results 

The  results  presented  so  far  for  small  k  and  v  show  that  the  variance  of  the  AVC 
estimator  is  weakly  distribution  dependent  and  tedious  to  compute.  We  derive 
asymptotic  variance  expressions  for  each  of  the  estimators  that  are  easy  to  estimate 
for  all  values  of  k. 
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3.5.1  Preliminaries 


Using  our  previous  notation  we  have 


pj  =  Pr{Xji  >  xu, w  /  j}. 


Define 


Yj(v)  =  number  of  wins  for  system  j  under  BEM  for  fixed  v 
Zj(v )  =  number  of  wins  for  system  j  under  AVC  for  fixed  v 


which  gives  us  point  estimators 


Pj(v) 
Pj(v ) 


Yj(v) 

v 

Zj(v)' 

vk 


Notice  these  are  the  same  estimators  we  defined  in  (3.1)  and  (3.4)  respectively, 
except  that  we  have  introduced  a  dependence  on  v.  Our  approach  is  based  on  the 
fact  that  standardized  versions  of  p  and  p  are  asymptotically  multivariate  normal 
(MVN). 


3.5.2  MLEs  and  AVC  Estimators 

Consider  the  asymptotic  behavior  of  the  MLEs.  It  is  well  known  that  as  v  — >  oo 

Vv(Pj  ~  Pj)  N  (0,pj(l  -  Pj)) .  (3.12) 

Similarly,  consider  the  asymptotic  behavior  of  our  AVC  estimator  as  v  goes  to 
infinity.  Following  Randles  and  Wolfe  (1979,  pp.  105-106)  let  i  be  an  integer  such 
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that  1  <  i  <  k  and  define 

H<?>  =  h<‘>(XIa, XtaJ 

and 

Hi?  =  h<‘>( 

where 

I  ajx  =  pj i  j  =  i 
\  ajl  ~f~  Pjl  j  7^  *• 

Then  define  the  covariance  terms 

d::“h,o,...,o  =  cov[4“>,4“»; 

=  E  [H\?H$]  -  pi 

where  the  only  1  in  the  subscript  of  sn  ui  i.u.  tj  is  in  the  ith  position  and  a  6 
{1  With  this  notation  a  represents  the  population  whose  parameter  we 

are  estimating,  and  i  represents  the  only  population  with  a  common  observation 
in  and  H$ .  Then  using  our  kernel  from  (3.5)  we  have  for  a  =  i 

=  Pr  {A'a„.  >  max{A-to„  A'(„,}} 

and  for  a  ^  i 

Ef/f(a) /f(°)]  =  Pr  1  ^ a<Xa  ^  ^ia‘  ’  ^aa<l  ^  }j  1 

41  i2  ~  \  Xapa>Xiai,Xapa>ma Xi^a,i{Xept}  ) 

This  leaves  us  with  just  two  covariance  expressions  for  each  population 

£<“>(l)  =  Co v[H%\H<g]  (3.13) 

and 

£W(2)  =  Cov[//<1“>,4“)],  a  +  i.  (3.14) 
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Let  N  —  kv.  We  then  state  the  following  theorem  due  to  Lehmann  (1963).  Some 
of  the  notation  has  been  simplified  for  our  context. 


Theorem  3.1  (Asymptotic  Distribution  of  AVC  Estimators) 

Let  Ua(X n, . . . ,  Alt); . . . ;  AAi, . . . ,  X^v)  =  pa  be  a  k-sample  U-statistic  for  the 
parameter  pa  of  degree  (1,...,1).  If  lim„_Kxl(?;/./V)  =  0  <  A*  <  1,  for  i  — 

1, . . . ,  k,  and  if  . . . ,  A*i)}2]  <  oo,  then  \/N(pa  —  pa)  has  a  limiting 

normal  distribution  with  mean  0  and  variance 


provided  a2  >  0. 

Proof:  See  Lehmann  (1963). 

Using  (3.13)  and  (3.14)  we  can  simplify  (3.15)  as 

<72  =  /t[Cw(l)  +  (/t-l)£w(2)]. 


(3,15) 


(3.16) 


From  Theorem  3.1  we  have  as  v  — >  oo 

V N(Pj  -  Pj)  ^  N(0,  a2)  j  =  1, . . . , k  (3.17) 

where  N  =  kv  and  a2  is  defined  in  (3.16).  In  comparing  a2  with  the  asymptotic 
MLE  variance  in  (3.12),  we  notice  that  the  MLE  expression  is  defined  as  yjv  — >  oo 
while  the  AVC  expression  is  defined  as  Vkv  — >  oo.  To  allow  a  fair  comparison 
between  these  asymptotic  variances,  we  multiply  (3.12)  through  by  \Jk,  to  obtain 

Vk(Vv(Pj  -pj))  N(0,  kpj(l  -pj))  (3.18) 
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as  y/kv  — >  oo.  We  now  have  comparable  asymptotic  variances  for  AVC  and 
MLE  from  (3.17)  and  (3.18),  respectively.  In  order  to  compare  these  asymptotic 
variances  with  the  exact  variance  of  pk  for  any  k  we  divide  each  expression  by  N 


to  give  us  the  following 

Var  ,4  (pj)  = 

\pA1  -pA 

(3.19) 

Var  A(pj)  w 

-[£(%)  +  (Jfc  -  1)£(J)(2)1. 

V 

(3.20) 

To  compare  these  results  we  remove  the  dependence  of  the  pjS  on  v  by  selecting 
a  value  for  9.  Figure  3.4  plots  v  ■  Var^(p2),  v  ■  Var ^(p2)  and  v  times  the  exact 
variance  found  in  (3.10)  against  v  for  k  =  2  exponential  populations  with  9  =  1.2. 
This  figure  illustrates  some  very  important  facts  about  the  relationship  among  the 
variances  of  our  estimators.  First,  notice  the  significant  reduction  in  the  asymptotic 
approximation  for  our  AVC  variance  over  the  MLE  variance.  Also  notice  how 
quickly  the  exact  AVC  variance  approaches  the  AVC  asymptotic  approximation. 
At  v  =  15,  the  difference  is  only  about  3%.  This  indicates  that  the  asymptotic 
approximation  for  the  AVC  variance  is  quite  good  at  relatively  small  values  of  v 
( v  >  15)  for  9  =  1.2.  Calculations  for  larger  values  of  9  at  v  =  15  show  the 
difference  between  the  exact  and  approximate  AVC  variance  still  about  3%  for 
9  =  2.0  and  increasing  to  a  difference  of  roughly  4%  at  9  =  3.0.  This  shows 
that  the  Var ^(p2)  provides  a  better  approximation  to  Var(p2)  when  the  difference 
between  systems  is  smaller. 

With  the  favorable  results  presented  above  regarding  the  accuracy  of  our 
asymptotic  approximation  for  AVC  variance,  we  incorporate  this  approximate 
variance  to  calculate  confidence  intervals  about  our  pjS.  Because  of  the  ease  in 
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1  3  5  7  9  11  13  15 


v 

Figure  3.4:  t;-Var(p2),  u-Var (p2)  and  f-Var^(p2)  for  k  =  2  Exponential  Populations 
with  6  =  1.2 

calculating  the  exact  variance  for  the  MLE  estimators,  we  use  this  exact  variance 
in  constructing  confidence  intervals  for  the  pjS. 

Suppose  we  have  generated  v  independent  observations  from  each  of  k  popu¬ 
lations.  We  compute  pj  and  pj  (j  =  1  using  (3.1)  and  (3.4),  respectively. 

We  then  estimate  Var ^ipj)  an(l  Var a(Pj)  using  (3.19)  and  (3.20),  respectively.  We 
perform  M  macro-replications  of  these  kv  observations  and  let  pJU  pji,  Var Aipj)i 
and  Var A{Pj)i  represent  the  point  estimates  and  variances  from  the  ith  macro¬ 
replication.  The  sample  variance  is  an  unbiased  estimator  of  the  true  variance  so 
we  calculate 
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Var  s(Pj) 


Eth  ( pj  -  pji} 
M  -  1 


where 


M  p.. 

-p  =  y ^ 

i=i 

We  can  also  calculate  average  variances  defined  as 


V=4(ft)  =  E VarA(Pi>< 

i= 1 


M 


vs^(ft)  =  EVaMfe)i 

i=l 


M 


We  then  estimate  the  bias  of  our  AVC  estimator  variance  as 


Bias 


=  Var A(pj)  -  Var s(pj). 


Let  Za/2  denote  the  a/2  quantile  of  the  standard  normal  distribution.  In 
constructing  our  confidence  intervals,  we  use  the  appropriate  variance  calculated 
in  each  macro-replication  and  then  using  normal  approximations  compute  for  the 


MLE 

Pji  ±  Za/2yJVaxA{pj)i 

and  for  AVC 

Pji  ±  Zai2\IVaxA(pj)i. 


We  then  count  how  many  of  the  M  intervals  formed  with  each  method  capture  pj 
and  compute  the  average  confidence  interval  width. 

We  provide  empirical  results  comparing  the  intervals  between  the  two  methods 
in  the  following  section. 
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3.6  Empirical  Study 


For  our  empirical  study  we  focus  on  the  difference  in  the  variance  of  the  MLEs 
as  compared  to  the  variance  of  our  AVC  estimators.  These  comparisons  take  two 
different  forms.  First  we  simply  examine  the  differences  in  the  variances  between 
the  two  methods.  Since  our  AVC  variance  estimator  is  an  approximation,  we  also 
obtain  an  unbiased  estimator  of  the  true  AVC  variance  and  compute  an  estimated 
bias  for  the  AVC  variance  estimator.  The  variance  reduction  obtained  is  shown  as 
Var5(pfc)/Var^(p/t).  We  then  construct  confidence  intervals  for  the  MLEs  and  AVC 
estimates  and  compare  the  percentage  of  coverage  and  average  confidence  interval 
width.  All  results  are  computed  for  pk  =  p^  using  exponential  population  distri¬ 
butions  as  described  in  §3.4.2.  Results  are  based  on  M=10,000  macro-replications 
for  each  of  the  v  values  listed. 

Table  3.1  presents  the  variance  results  for  each  method  for  k  =  2  to  5  popula¬ 
tions  and  a  number  of  different  values  of  v.  We  start  off  with  v  =  20  to  see  how 
well  our  approximation  works  for  small  values  of  v.  We  find  a  significant  negative 
bias  in  our  AVC  estimator  variance  at  v  =  20  for  all  V  At  v  =  50  this  bias  drops 
to  between  6%  and  7%  for  k  =  2  and  less  than  10%  for  k  =  3.  At  k  =  5  the 
bias  is  stll  more  than  10%.  Moving  up  to  v  =  100,  the  bias  for  both  k  —  2  and 
3  drops  below  7%  and  for  k  =  5  the  bias  drops  below  10%.  At  v  =  200  the  bias 
effectively  goes  to  zero  for  all  k.  These  results  indicate  that  we  can  substantially 
underestimate  Var (pk)  by  using  Var^(pfe)  for  v  <  50.  The  last  column  of  Table  3.1 
shows  the  ratio  of  Vars^/^/Var^p*,)  to  illustrate  the  variance  reduction  with  the 
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AVC  estimator.  Focusing  on  the  results  for  v  =  200  at  9*  —  1.2,  we  see  roughly 
a  33%  reduction  at  k  =  2,  roughly  a  42%  reduction  at  k  =  3,  and  roughly  a  45% 
reduction  at  k  =  5. 

Looking  at  the  95%  confidence  interval  results  in  Table  3.2,  we  see  poor  coverage 
and  relatively  large  intervals  for  v  <  50  for  both  MLE  and  AVC.  This  indicates 
that  the  normal  approximation  is  not  particularly  good  for  either  method  at  small 
values  of  v.  At  v  =  100  both  MLE  and  AVC  coverage  jump  up  to  94%  to  95% 
in  nearly  all  cases,  with  MLE  coverage  slightly  better  than  the  AVC  coverage. 
However,  at  v  —  200  the  AVC  coverage  slightly  exceeds  the  MLE  coverage,  with 
both  methods  achieving  95%  coverage  in  almost  all  case.  In  all  cases  the  average 
AVC  interval  width  shows  a  20-25%  reduction  over  the  average  MLE  interval  width. 

3.7  Conclusions 

We  have  shown  that  our  AVC  estimator  for  multinomial  success  probabilities  has 
a  no  larger  variance  than  the  MLE.  Our  results  indicate  a  reduction  in  variance  on 
the  order  of  33%  for  k  =  2,  increasing  to  roughly  45%  for  k  =  5  populations  can  be 
expected.  We  have  also  shown  that  this  variance  reduction  depends  weakly  on  the 
distribution  of  the  performance  measures  from  the  populations.  An  asymptotic 
approximation  for  our  AVC  estimator  variance  was  found  and  this  was  used  to 
compare  to  the  MLE  variance  directly,  and  through  the  construction  of  confidence 
intervals.  For  small  v  we  saw  that  our  AVC  estimator  variance  has  a  relatively 
large  negative  bias.  However,  for  v  >  100,  this  bias  became  insignificant.  The 
confidence  interval  results  show  similar  percent  coverage  with  AVC  and  MLE  for 
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k 

6* 

V 

Var A(pk) 

Bias(Var^(pfc)) 

Var  s(Pk) 
Var  (vie) 

1.2 

20 

.0086 

.0074 

.0118 

-.0012 

.7288 

2 

1.2 

50 

.0034 

.0032 

.0049 

-.0002 

.6939 

1.2 

100 

.0017 

.0016 

.0025 

-.0001 

.6800 

1.2 

200 

.0008 

.0008 

.0012 

.0000 

.6667 

.0076 

.0105 

-.0011 

.7238 

2 

50 

.0030 

.0028 

.0044 

-.0002 

.6818 

.0015 

.0014 

.0022 

-.0001 

.6818 

£ 

.0007 

.0007 

.0011 

.0000 

.6364 

1.2 

■ 

-.0014 

.6757 

3 

1.2 

1 

I 

-.0003 

.6522 

1.2 

.0014 

-.0001 

.6522 

1.2 

.0007 

.0000 

.5833 

2.0 

20 

.0085 

.0069 

.0119 

-.0016 

.7143 

3 

2.0 

50 

.0034 

.0031 

.0049 

-.0003 

.6939 

2.0 

100 

.0017 

.0016 

.0025 

-.0001 

.6800 

2.0 

200 

.0008 

.0008 

.0012 

.0000 

.6667 

1.2 

20 

.0053 

.0038 

.0084 

-.0015 

.6310 

5 

1.2 

50 

.0021 

.0018 

.0035 

-.0003 

.6000 

1.2 

100 

.0011 

.0010 

.0018 

-.0001 

.6111 

1.2 

200 

.0005 

.0005 

.0009 

.0000 

.5555 

2.0 

20 

.0073 

.0053 

.0105 

-.0020 

.6952 

5 

2.0 

50 

.0029 

.0026 

.0043 

-.0003 

.6744 

2.0 

100 

.0015 

.0014 

.0022 

-.0001 

.6818 

2.0 

200 

.0007 

.0007 

.0011 

.0000 

.6364 

Table  3.1:  AVC  and  MLE  Variance  for  pt 
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AVC 

MLE 

k 

e* 

V 

Percentage 

Coverage 

Average 

Width 

Percentage 

Coverage 

Average 

Width 

1.2 

20 

.92 

.34 

.92 

.42 

2 

1.2 

50 

.94 

.22 

.93 

.27 

1.2 

100 

.94 

.16 

.94 

.19 

1.2 

200 

.95 

.11 

.94 

.14 

.91 

.31 

.92 

.40 

2 

2.0 

.94 

.21 

.95 

.26 

2.0 

.94 

.15 

.96 

.18 

2.0 

.95 

.10 

.95 

.13 

1.2 

.90 

.31 

.94 

.41 

3 

EH 

.93 

.20 

.94 

.27 

1.2 

.94 

.15 

.95 

.19 

EH 

.95 

.10 

.94 

.13 

2.0 

20 

.91 

.32 

.96 

.43 

3 

2.0 

50 

.93 

.22 

.93 

.27 

2.0 

100 

.94 

.16 

.94 

.20 

2.0 

200 

.95 

.11 

.95 

.14 

1.2 

20 

.87 

.24 

.94 

.35 

5 

1.2 

50 

.92 

.17 

.95 

.23 

1.2 

100 

.93 

.12 

.95 

.16 

1.2 

200 

.94 

.09 

.94 

.12 

2.0 

20 

.88 

.28 

.93 

.40 

5 

2.0 

50 

.92 

.20 

.94 

.26 

2.0 

100 

.94 

.14 

.95 

.18 

2.0 

200 

.95 

.10 

.95 

.13 

Table  3.2:  AVC  and  MLE  95%  Confidence  Intervals  for 
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v  >  100  with  the  average  AVC  interval  width  20-25%  less  than  MLE.  Our  AVC 
estimators  make  more  efficient  use  of  the  data  already  available  to  provide  a  more 
precise  set  of  estimators  for  the  multinomial  success  probabilities. 
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Appendix  A 


PROOFS 


Proof  of  Theorem  2.1: 

From  (2.9),  as  v  — >  oo  we  have  Pi(v)  1  jk  for  all  i.  Then  from  Lemma  2.1 
1  A^{v)  N 

V  Mv)  ) 


(  ~S  ) 

f  l/k(l  -  1/k) 

-1/fc2 

- 1/k 2  \ 

" 

MVN 

9 

-1/k2 

l/k(l-l/k)  •• 

-1/A;2 

{-*) 

{  - 1/k 2 

—  1/k2 

••  iA(i-i/fc)J 

Since  Ai(v )  —  Ai(v)  =  Ai(v )  +  X^_2  A(w)  (®  ^  1)>  we  take  the  difference  of  MVN 
random  variables  and  obtain 

f  A2(f)  -  Ai(r>)  > 

V  -  Ai(v)  J 


-M  1 

/  1 

1/2  • 

1/2  \ 

MVN 

,2/A; 

1/2 

1  • 

'•  1/2 

V  1/2 

1/2  •• 

•  1  ) 

(A.l) 


78 


Assuming  population  1  is  the  best,  in  terms  of  PCS  we  can  state 


PCS  =  Urn_  Pr  {!',(>>)  >  1»,  Vj  A  1} 

=  lim  pr  f YM-v/k  >  YM-vJk 

=  lim  Pr{Ai(u)  >  Aj(v),  Vj  #  1} 


,  Vj#  1 


=  lim  Pr{Aj(n)  —  A^u)  <  0,  V/'  #  1} 
=  Pr{VPj  <  0,  j  =  2,3, 


(A. 2) 


where  (W2, . . . ,  M4)  ~  (A.l).  If  we  add  M  to  each  Wj  to  obtain  a  random  vector 
with  a  mean  of  zero,  then  from  (A. 2)  we  have 


pcT  -  Pr m±M<  *  ,-_2 . * 

l  AA  AA 


=  Pr  r  Qj  <  —^==,  j  =  2,...,k 


where 


MVN 


Aa 


1  1/2  •••  1/2 

1/2  1  •••  1/2 


1/2  1/2 


Proof  of  Theorem  2.2: 

Define  X(IV),  the  covariance  matrix  computed  for  {F\v\  . . . ,  F^),  as 


£(7V)  =  Var 


Pi  (AO 

Pk(N) 


Let  E N[Pj{N)]  =  Pj{N)  where  Ejv[Pj(A0]  denotes  the  expected  value  at  sample 
size  N  (N  =  kv).  Lehmann  (1963,  pp.  964-965)  shows  that  Lemma  2.2  holds  even 


79 


if  the  distributions  of  the  data  depend  upon  the  sample  size  provided 

£(AQ  — ►  E, 

where  E  is  nonsingular  as  N  — >  oo.  Since  =  1,  we  work  with  p2, ...  ,pk 

to  be  nonsingular.  We  assume 

/  y/N(p2(N)-p2(N))  \ 

Var  :  — » E 

V  y/N(pk(N)-pk(N)) 

for  any  reasonable  set  of  F-v> .  We  define  E  as  ||cda,6)j|  in  (2.12). 

To  proceed,  we  need  to  consider  the  covariance  terms  defined  in  Theorem  2.2. 

?<“’*>  =  Elflft'ifJJ’]  - 


since 

rr(a)  =  (  1,  if*  aaQ  >  Xfo,  W  ^  a 

41  —  (  0,  otherwise 

rr(6)  _  f  1,  ifXbf}b>XepcV£^b 

12  (0,  otherwise 


Using  our  kernel  from  (2.11)  we  can  express 

Elflfftf&’l  =  E[  n  4>(Xaa.  -  Xm)  x  n  </>(X.c„  -  Xm)l 

We  can  then  write  this  as  a  probability  statement  combining  the  two  indices  to 
come  up  with 

=  Pr  {a„„„  >  max{A',„„  A'w,}}  . 
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This  particular  case  is  easy  to  illustrate  and  we  extend  this  development  for  other 
expected  value  terms  in  less  detail. 

There  are  a  number  of  different  cases  we  need  to  consider  for  the  covariance 
terms.  In  our  notation  £fa,b\  a  represents  the  population  with  the  largest  value  in 
H\a^ ;  b  represents  the  population  with  the  largest  value  in  H^;  and  i  represents  the 
one  population  that  has  the  same  observation  in  both  vectors  (pseudo-replications). 
We  can  enumerate  the  different  cases  for  the  covariance  terms  based  on  the  values 
of  a,  b  and  i.  We  have  the  following  four  cases. 


1.  a  —  b  =  i\ 

2-  a  =  b±  *;  d°’a) 

3.  ajLb^i;  ^a,b) 

4.  a  b,  a  =  %  or  b  =  i;  ^a,6l  or 


Asymptotically  Xj  ~  F  Vj  which  allows  us  to  construct  distribution  free  ex¬ 
pressions  for  £}a,b'>  since  each  Xtj  has  the  same  probability  of  being  the  largest  value 
in  a  single  vector.  For  case  1  we  have 


&a)  =  Pr  {Xa«a  >  max{Xtef,  AT^f}|  -  ~ 

1  1 
2k  -  1  _ 


(A.3) 


This  follows  since  we  have  2(k  —  1)  +  1  independent  and  identically  distributed 
random  variables  and  we  want  the  probability  that  a  particular  one  is  the  largest. 
In  our  context  this  means  we  want  the  probability  that  a  single  observation  from 
population  a  is  a  winner  in  two  separate  vectors  containing  no  other  common 
observations. 


81 


Case  2  is  more  difficult  to  approach.  We  have 


t(a-a)  _  pr  |  j^«a„  >  X%ai  i  Xaaa  >  max^a)j{A^af}; 

Xa/3a  A ia;  5  A a/3a  ^ 

where  the  common  observation  in  each  of  the  two  vectors  is  not  the  largest,  and 
both  observations  from  population  a  are  the  largest  in  their  respective  vectors. 
The  ordering  of  the  random  variables  from  this  pair  of  vectors  must  look  like  the 
following: 

Xaaa  >  {ho  f  the  X(Qt }  >  Xa0a  > 

{{k  —  2)  -  h of  the Xta( ,  (fc  —  2)  of  the X^n  and  Xicei} 

or  interchanging  X0Qa  and  Xa/3a 

Xada  >  {h  of  the  Xep( }  > 

^ acta  ^ 

{{k  —  2)  —  h of  the  X^n  (A:  —  2)  of  the  A&*p  and  Xiai} 
where  h  =  0, 1,. . .  ,k  —  2.  For  each  subset  of  size  h  there  are 

h\{(k  —  2)  —  h  +  (k  —  2)  +  1)\ 


equally  likely  orderings  and  there  are 

k  —  2 
h 

ways  to  select  an  h.  With  a  total  of  (2 k  —  1)!  possible  orderings  we  then  have 

2Ej;2„(*T2)«((fc-2)-/>+(fc-2)  +  l)!  j 

^  ’  =  {2k -1)1  k 2 

2  1 

k(2k  —  1)  k2'  ' 
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=  k  Y.C 


{a, a) 


i=  1 


=  k  [?(“■“> + (k  ~  i)^ 

1  1 


k 


2k  -  1  k? 
k-  1 


+  {k-  1) 


k(2k  —  1)  A;2 


(A.7) 


2k  -  1 

With  the  covariance  expressions  from  (A. 5)  and  (A. 6)  we  can  find  the  off- 
diagonal  terms  of  from  (2.12)  as 


a 


(a,b) 


£{a,b)  _j_  £(a,6)  ^  g( a’b ) 

£^a,b 


k 


2  l 

-1 


1 


A;(2A;  -  1)  k°- 


~  T? 1  +  (k  -  2) 


k(2k  —  1)  A;2/J 


(A.8) 


2k  -  1 

Combining  the  terms  in  (A.8)  with  those  in  (A.7)  we  let  rj  =  1  /(k  —  1)  and  we 
have 


/  1 


£  = 


'  k-  r 

2k  -  1 


■v  ■■■  -v\ 


\  -V  -V 

From  Lehman  (1963)  we  know  that 


??  1  •••  -rj 

•  1  / 


(A. 9) 


\ 

(0) 

==>  MVN 

/ 

\0) 

(A. 10) 


Asymptotically  each  pj  =  l/k  with  our  model,  so  each  difference  in  (A. 10)  becomes 
y/N(pj(N)  —  1/A;).  The  resultant  shift  in  the  mean  of  (A. 10)  can  be  found  as 

VN{pj(N)  -  l/k)  =  VNipjW-il/k-S/y/dft-VNS/yfi 
=  VN(pj{N)  -  (l/k)  -  6/yfi))  -  VkS 
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Then 


PCS 


=  Pr 


kVkS 


I  Wj  +  ky/kS  ^ 

\  \j2k/ (2k  —  1)  y/2k/(2k  -  1) 


i  =  2, 


Pr 


kS 

^2/(2*  -  1)  ’ 


where 


□ 


'l 

r°'l 

~  MVN 

\  Qk  ) 

/  1 

1/2 

•  1/2  \ 

1/2 

1  •• 

•  1/2 

V  1/2 

1/2 

•  1  ) 

Appendix  B 


TABLES  OF  RESULTS 


The  following  tables  are  based  on  100,000  macro-replications.  The  ratio  column 
is  PCSavc/PCSbem  and  the  s.e.  columns  are  the  respective  standard  errors.  The 
population  distributions  used  are  discussed  in  §2.6. 


V 

PCSbem 

PCSavc 

Ratio 

s.e.  PCSbem 

s.e.  PCSavc 

2 

0.5463 

0.5566 

1.0188 

0.0016 

0.0016 

4 

0.5684 

0.5849 

1.0290 

0.0016 

0.0016 

6 

0.5853 

0.6058 

1.0350 

0.0016 

0.0015 

8 

0.5998 

0.6248 

1.0416 

0.0015 

0.0015 

10 

0.6106 

0.6370 

1.0432 

0.0015 

0.0015 

12 

0.6218 

0.6502 

1.0457 

0.0015 

0.0015 

14 

0.6317 

0.6615 

1.0471 

0.0015 

0.0015 

16 

0.6405 

0.6722 

1.0495 

0.0015 

0.0015 

18 

0.6475 

0.6813 

1.0523 

0.0015 

0.0015 

20 

0.6559 

0.6901 

1.0521 

0.0015 

0.0015 

22 

0.6641 

0.6981 

1.0512 

0.0015 

0.0015 

24 

0.6711 

0.7057 

1.0515 

0.0015 

0.0014 

26 

0.6770 

0.7133 

1.0536 

0.0015 

0.0014 

28 

0.6825 

0.7215 

1.0572 

0.0015 

0.0014 

30 

0.6893 

0.7277 

1.0558 

0.0015 

0.0014 

32 

0.6944 

0.7331 

1.0557 

0.0015 

0.0014 

34 

0.7004 

0.7386 

1.0546 

0.0014 

0.0014 

36 

0.7053 

0.7454 

1.0567 

0.0014 

0.0014 

38 

0.7110 

0.7520 

1.0577 

0.0014 

0.0014 

40 

0.7164 

0.7574 

1.0573 

0.0014 

0.0014 

42 

0.7210 

0.7635 

1.0590 

0.0014 

0.0013 

44 

0.7257 

0.7694 

1.0602 

0.0014 

0.0013 

46 

0.7299 

0.7727 

1.0586 

0.0014 

0.0013 

48 

0.7342 

0.7781 

1.0598 

0.0014 

0.0013 

50 

0.7378 

0.7832 

1.0615 

0.0014 

0.0013 

Table  B.l:  Simulation  Results  for  Exponential  Populations  with  9  =  1.2  and  k  =  2 
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PCSbem 

0.2853 

0.3032 

0.3150 

0.3258 

0.3335 

0.3427 

0.3492 

0.3573 

0.3640 

0.3716 

0.3795 

0.3839 

0.3893 

0.3949 

0.4000 

0.4050 

0.4093 

0.4147 

0.4184 

0.4245 

0.4287 

0.4327 

0.4374 

0.4421 

0.4465 


PCSavc 
0.298 
0.317 
0.335 
0.350 
0.361 
0.371 
0.381 
0.392 
0.401 
0.410 
0.420 
0.426 
0.435 
0.443 
0.451 
0.456 
0.463 
0.468 
0.475 
0.482 
0.488 
0.493 
0.499 
0.503 
0.51C 


Ratio 

1.0456 

1.0465 

1.0643 

1.0745 

1.0834 

1.0831 

1.0935 

1.0978 

1.1016 

1.1041 

1.1075 

1.1113 

1.1173 

1.1220 

1.1282 

1.1281 

1.1329 

1.1286 

1.1352 

1.1354 

1.1398 

1.1408 

1.1417 

1.1393 

1.1423 


s.e.  PCS 


0.0014 

0.0015 

0.0015 

0.0015 

0.0015 

0.0015 

0.0015 

0.0015 

0.0015 

0.0015 

0.0015 

0.0015 

0.0015 

0.0015 

0.0015 

0.0016 

0.0016 

0.0016 

0.0016 

0.0016 

0.0016 

0.0016 

0.0016 

0.0016 

0.0016 


s.e.  PCS 


Table  B.3:  Simulation  Results  for  Exponential  Populations  with  6  =  1.2  and  k 


V 

PCSbem 

PCSavc 

Ratio 

s.e.  PCSbem 

s.e.  PCSavc 

2 

0.2299 

0.2413 

1.0495 

0.0013 

0.0014 

4 

0.2443 

0.2581 

1.0566 

0.0014 

0.0014 

6 

0.2543 

0.2709 

1.0656 

0.0014 

0.0014 

8 

0.2637 

0.2849 

1.0806 

0.0014 

0.0014 

10 

0.2702 

0.2952 

1.0925 

0.0014 

0.0014 

12 

0.2786 

0.3047 

1.0938 

0.0014 

0.0015 

14 

0.2840 

0.3122 

1.0989 

0.0014 

0.0015 

16 

0.2898 

0.3217 

1.1098 

0.0014 

0.0015 

18 

0.2955 

0.3289 

1.1128 

0.0014 

0.0015 

20 

0.3029 

0.3371 

1.1129 

0.0015 

0.0015 

22 

0.3079 

0.3470 

1.1268 

0.0015 

0.0015 

24 

0.3130 

0.3524 

1.1261 

0.0015 

0.0015 

26 

0.3172 

0.3596 

1.1337 

0.0015 

0.0015 

28 

0.3224 

0.3669 

1.1382 

0.0015 

0.0015 

30 

0.3267 

0.3736 

1.1435 

0.0015 

0.0015 

32 

0.3305 

0.3787 

1.1458 

0.0015 

0.0015 

34 

0.3345 

0.3854 

1.1520 

0.0015 

0.0015 

36 

0.3390 

0.3901 

1.1507 

0.0015 

0.0015 

38 

0.3431 

0.3955 

1.1527 

0.0015 

0.0015 

40 

0.3477 

0.4018 

1.1558 

0.0015 

0.0016 

42 

0.3504 

0.4074 

1.1625 

0.0015 

0.0016 

44 

0.3555 

0.4128 

1.1612 

0.0015 

0.0016 

46 

0.3593 

0.4180 

1.1634 

0.0015 

0.0016 

48 

0.3632 

0.4220 

1.1619 

0.0015 

0.0016 

50 

0.3674 

0.4267 

1.1614 

0.0015 

0.0016 

Table  B.4:  Simulation  Results  for  Exponential  Populations  with  9  =  1.2  and  k  =  5 
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V 

PCSbem 

PCSavc 

Ratio 

s.e.  PCSbem 

s.e.  PCSavc 

2 

0.6675 

0.7003 

1.0491 

0.0015 

0.0014 

4 

0.7408 

0.7901 

1.0665 

0.0014 

0.0013 

6 

0.7903 

0.8428 

1.0665 

0.0013 

0.0012 

8 

0.8270 

0.8772 

1.0606 

0.0012 

0.0010 

10 

0.8546 

0.9040 

1.0578 

0.0011 

0.0009 

12 

0.8779 

0.9243 

1.0528 

0.0010 

0.0008 

14 

0.8964 

0.9387 

1.0472 

0.0010 

0.0008 

16 

0.9117 

0.9516 

1.0437 

0.0009 

0.0007 

18 

0.9239 

0.9604 

1.0395 

0.0008 

0.0006 

20 

0.9347 

0.9677 

1.0353 

0.0008 

0.0006 

22 

0.9439 

0.9739 

1.0318 

0.0007 

0.0005 

24 

0.9522 

0.9789 

1.0281 

0.0007 

0.0005 

26 

0.9577 

0.9826 

1.0260 

0.0006 

0.0004 

28 

0.9635 

0.9855 

1.0228 

0.0006 

0.0004 

30 

0.9680 

0.9882 

1.0208 

0.0006 

0.0003 

32 

0.9727 

0.9906 

1.0184 

0.0005 

0.0003 

34 

0.9762 

0.9923 

1.0165 

0.0005 

0.0003 

36 

0.9792 

0.9930 

1.0141 

0.0005 

0.0003 

38 

0.9816 

0.9943 

1.0130 

0.0004 

0.0002 

40 

0.9837 

0.9953 

1.0118 

0.0004 

0.0002 

42 

0.9856 

0.9959 

1.0104 

0.0004 

0.0002 

44 

0.9875 

0.9963 

1.0089 

0.0004 

0.0002 

46 

0.9891 

0.9973 

1.0083 

0.0003 

0.0002 

48 

0.9905 

0.9976 

1.0071 

0.0003 

0.0002 

50 

0.9914 

0.9978 

1.0064 

0.0003 

0.0001 

Table  B.6:  Simulation  Results  for  Exponential  Populations  with  9  =  2.0  and  k  =  2 
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V 

PCSbem 

PCSavc 

Ratio 

s.e.  PCSbem 

s.e.  PCSavc 

2 

0.5005 

0.5560 

1.1109 

0.0016 

0.0016 

4 

0.5936 

0.6481 

1.0917 

0.0016 

0.0015 

6 

0.6460 

0.7133 

1.1043 

0.0015 

0.0014 

8 

0.6868 

0.7635 

1.1116 

0.0015 

0.0013 

10 

0.7266 

0.8022 

1.1041 

0.0014 

0.0013 

12 

0.7579 

0.8345 

1.1011 

0.0014 

0.0012 

14 

0.7827 

0.8600 

1.0988 

0.0013 

0.0011 

16 

0.8082 

0.8823 

1.0917 

0.0012 

0.0010 

18 

0.8280 

0.9007 

1.0878 

0.0012 

0.0009 

20 

0.8457 

0.9157 

1.0827 

0.0011 

0.0009 

22 

0.8638 

0.9289 

1.0753 

0.0011 

0.0008 

24 

0.8756 

0.9376 

1.0707 

0.0010 

0.0008 

26 

0.8879 

0.9468 

1.0663 

0.0010 

0.0007 

28 

0.9003 

0.9554 

1.0612 

0.0009 

0.0007 

30 

0.9097 

0.9609 

1.0562 

0.0009 

0.0006 

32 

0.9185 

0.9666 

1.0525 

0.0009 

0.0006 

34 

0.9260 

0.9718 

1.0495 

0.0008 

0.0005 

36 

0.9334 

0.9757 

1.0453 

0.0008 

0.0005 

38 

0.9392 

0.9790 

1.0424 

0.0008 

0.0005 

40 

0.9444 

0.9822 

1.0400 

0.0007 

0.0004 

42 

0.9507 

0.9842 

1.0352 

0.0007 

0.0004 

44 

0.9547 

0.9866 

1.0334 

0.0007 

0.0004 

46 

0.9590 

0.9887 

1.0310 

0.0006 

0.0003 

48 

0.9629 

0.9901 

1.0282 

0.0006 

0.0003 

50 

_ 

0.9670 

0.9908 

1.0247 

0.0006 

0.0003 

Table  B.7:  Simulation  Results  for  Exponential  Populations  with  9  =  2.0  and  k  =  3 


Table  B.8:  Simulation  Results  for  Exponential  Populations  with  0  =  2.0  and  k  =  4 
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V 

PCSbem 

PCSavc 

Ratio 

s.e.  PCSbem 

s.e.  PCSavc 

2 

0.3840 

0.0015 

0.0015 

4 

0.4069 

0.4634 

1.1388 

0.0016 

0.0016 

6 

0.4509 

0.5271 

1.1689 

0.0016 

0.0016 

8 

0.4941 

0.5816 

1.1771 

0.0016 

0.0016 

10 

0.5307 

0.6256 

1.1789 

0.0016 

0.0015 

12 

0.5635 

0.6674 

1.1843 

0.0016 

0.0015 

14 

0.5931 

0.7017 

1.1831 

0.0016 

0.0014 

16 

0.6210 

0.7338 

1.1817 

0.0015 

0.0014 

18 

0.6474 

0.7610 

1.1754 

0.0015 

0.0013 

20 

0.6708 

0.7857 

1.1713 

0.0015 

0.0013 

22 

0.6923 

0.8085 

1.1680 

0.0015 

0.0012 

24 

0.7113 

0.8273 

1.1631 

0.0014 

0.0012 

26 

0.7300 

0.8431 

1.1548 

0.0014 

0.0012 

28 

0.7472 

0.8592 

1.1498 

0.0014 

0.0011 

30 

0.7633 

0.8720 

1.1424 

0.0013 

0.0011 

32 

0.7785 

0.8847 

1.1364 

0.0013 

0.0010 

34 

0.7905 

0.8952 

1.1325 

0.0013 

0.0010 

36 

0.8045 

0.9054 

1.1254 

0.0013 

0.0009 

38 

0.8158 

0.9148 

1.1214 

0.0012 

0.0009 

40 

0.8274 

0.9230 

1.1156 

0.0012 

0.0008 

42 

0.8379 

0.9314 

1.1116 

0.0012 

0.0008 

44 

0.8483 

0.9381 

1.1059 

0.0011 

0.0008 

46 

0.8562 

0.9441 

1.1027 

0.0011 

0.0007 

48 

0.8652 

0.9479 

1.0956 

0.0011 

0.0007 

50 

0.8732 

0.9537 

1.0922 

0.0011 

0.0007 

Table  B.9:  Simulation  Results  for  Exponential  Populations  with  9  =  2.0  and  k 
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V 

PCSbem 

PCSavc 

Ratio 

s.e.  PCSbem 

s.e.  PCSavc 

2 

0.5468 

0.5601 

1.0244 

0.0016 

0.0016 

4 

0.5688 

0.5881 

1.0338 

0.0016 

0.0016 

6 

0.5849 

0.6064 

1.0368 

0.0016 

0.0015 

8 

0.5993 

0.6221 

1.0381 

0.0016 

0.0015 

10 

0.6122 

0.6363 

1.0395 

0.0015 

0.0015 

12 

0.6220 

0.6490 

1.0433 

0.0015 

0.0015 

14 

0.6322 

0.6617 

1.0466 

0.0015 

0.0015 

16 

0.6400 

0.6704 

1.0474 

0.0015 

0.0015 

18 

0.6479 

0.6793 

1.0484 

0.0015 

0.0015 

20 

0.6585 

0.6900 

1.0478 

0.0015 

0.0015 

22 

0.6647 

0.6989 

1.0514 

0.0015 

0.0015 

24 

0.6715 

0.7073 

1.0533 

0.0015 

0.0014 

26 

0.6779 

0.7136 

1.0527 

0.0015 

0.0014 

28 

0.6844 

0.7225 

1.0557 

0.0015 

0.0014 

30 

0.6913 

0.7283 

1.0535 

0.0015 

0.0014 

32 

0.6953 

0.7349 

1.0570 

0.0015 

0.0014 

34 

0.7017 

0.7410 

1.0561 

0.0014 

0.0014 

36 

0.7066 

0.7485 

1.0593 

0.0014 

0.0014 

38 

0.7115 

0.7538 

1.0595 

0.0014 

0.0014 

40 

0.7167 

0.7572 

1.0565 

0.0014 

0.0014 

42 

0.7214 

0.7636 

1.0584 

0.0014 

0.0013 

44 

0.7252 

0.7679 

1.0588 

0.0014 

0.0013 

46 

0.7305 

0.7731 

1.0582 

0.0014 

0.0013 

48 

0.7338 

0.7787 

1.0612 

0.0014 

0.0013 

50 

0.7376 

0.7824 

1.0607 

0.0014 

0.0013 

Table  B.ll:  Simulation  Results  for  Uniform  Populations  with  9  =  1.2  and  k  =  2 
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V 

PCSbem 

PCSavc 

Ratio 

s.e.  PCSbem 

s.e.  PCSavc 

2 

0.3759 

0.4030 

1.0720 

0.0015 

0.0016 

4 

0.3992 

0.4197 

1.0514 

0.0015 

0.0016 

6 

0.4133 

0.4372 

1.0580 

0.0016 

0.0016 

8 

0.4236 

0.4534 

1.0704 

0.0016 

0.0016 

10 

0.4351 

0.4643 

1.0671 

0.0016 

0.0016 

12 

0.4459 

0.4795 

1.0753 

0.0016 

0.0016 

14 

0.4536 

0.4909 

1.0820 

0.0016 

0.0016 

16 

0.4644 

0.5007 

1.0783 

0.0016 

0.0016 

18 

0.4717 

0.5109 

1.0830 

0.0016 

0.0016 

20 

0.4793 

0.5215 

1.0880 

0.0016 

0.0016 

22 

0.4877 

0.5314 

1.0896 

0.0016 

0.0016 

24 

0.4922 

0.5406 

1.0983 

0.0016 

0.0016 

26 

0.4977 

0.5479 

1.1008 

0.0016 

0.0016 

28 

0.5061 

0.5570 

1.1006 

0.0016 

0.0016 

30 

0.5127 

0.5631 

1.0984 

0.0016 

0.0016 

32 

0.5177 

0.5709 

1.1028 

0.0016 

0.0016 

34 

0.5229 

0.5774 

1.1042 

0.0016 

0.0016 

36 

0.5288 

0.5857 

1.1077 

0.0016 

0.0016 

38 

0.5338 

0.5919 

1.1088 

0.0016 

0.0016 

40 

0.5397 

0.5975 

1.1071 

0.0016 

0.0016 

42 

0.5443 

0.6042 

1.1102 

0.0016 

0.0015 

44 

0.5487 

0.6097 

1.1111 

0.0016 

0.0015 

46 

0.5533 

0.6156 

1.1126 

0.0016 

0.0015 

48 

0.5590 

0.6205 

1.1099 

0.0016 

0.0015 

50 

0.5619 

0.6254 

1.1131 

0.0016 

0.0015 

Table  B.12:  Simulation  Results  for  Uniform  Populations  with  9  =  1.2  and  k 
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PCSbem 

0.2866 

0.3066 

0.3150 

0.3269 

0.3342 

0.3443 

0.3506 

0.3592 

0.3645 

0.3721 

0.3792 

0.3848 

0.3891 

0.3951 

0.4005 

0.4057 

0.4100 

0.4148 

0.4197 

0.4233 

0.4286 

0.4340 

0.4370 

0.4415 

0.4464 


PCSavc 

0.3141 

0.3237 

0.3416 

0.3527 

0.3658 

0.3793 

0.3889 

0.3963 

0.4066 

0.4155 

0.4252 

0.4320 

0.4394 

0.4462 

0.4550 

0.4604 

0.4675 

0.4743 

0.4799 

0.4859 

0.4927 

0.4964 

0.5037 

0.5075 

0.5119 


Ratio 


959 
1558 
1846 
1786 
1945 
017 
1.1092 
1.1034 
1.1154 
1.1166 
1.1212 
1.1228 
1.1295 
1.1291 
1.1362 
1.1348 
1.1402 
1.1433 
1.1435 
1.1477 
1.1494 
1.1439 
1.1527 
1.1496 
1.1467 


s.e.  PCSDem 


0. 

0. 

0. 

0. 

0. 

0. 
0.0015 
0.0015 
0.0015 
0.0015 
0.0015 
0.0015 
0.0015 
0.0015 
0.0016 
0.0016 
0.0016 
0.0016 
0.0016 
0.0016 
0.0016 
0.0016 
0.0016 
0.0016 
0.0016 


s.e.  PCS 


0.0015 

0.0015 

0.0016 

0.0016 

0.0016 

0.0016 

0.0016 

0.0016 

0.0016 

0.0016 

0.0016 

0.0016 

0.0016 

0.0016 

0.0016 

0.0016 

0.0016 

0.0016 

0.0016 


Table  B.13:  Simulation  Results  for  Uniform  Populations  with  9  =  1.2  and  k 
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V 

PCSbem 

PCSavc 

Ratio 

s.e.  PCSbem 

s.e.  PCSavc 

2 

0.2307 

0.2566 

1.1122 

0.0013 

0.0014 

4 

0.2459 

0.2660 

1.0815 

0.0014 

0.0014 

6 

0.2544 

0.2810 

1.1045 

0.0014 

0.0014 

8 

0.2631 

0.2914 

1.1074 

0.0014 

0.0014 

10 

0.2703 

0.3020 

1.1172 

0.0014 

0.0015 

12 

0.2782 

0.3130 

1.1250 

0.0014 

0.0015 

14 

0.2843 

0.3197 

1.1247 

0.0014 

0.0015 

16 

0.2922 

0.3278 

1.1218 

0.0014 

0.0015 

18 

0.2959 

0.3356 

1.1342 

0.0014 

0.0015 

20 

0.3029 

0.3437 

1.1345 

0.0015 

0.0015 

22 

0.3075 

0.3530 

1.1480 

0.0015 

0.0015 

24 

0.3129 

0.3597 

1.1495 

0.0015 

0.0015 

26 

0.3154 

0.3651 

1.1578 

0.0015 

0.0015 

28 

0.3224 

0.3714 

1.1522 

0.0015 

0.0015 

30 

0.3259 

0.3788 

1.1624 

0.0015 

0.0015 

32 

0.3298 

0.3833 

1.1621 

0.0015 

0.0015 

34 

0.3346 

0.3903 

1.1666 

0.0015 

0.0015 

36 

0.3405 

0.3960 

1.1632 

0.0015 

0.0015 

38 

0.3436 

0.4020 

1.1699 

0.0015 

0.0016 

40 

0.3463 

0.4071 

1.1753 

0.0015 

0.0016 

42 

0.3519 

0.4125 

1.1720 

0.0015 

0.0016 

44 

0.3564 

0.4161 

1.1673 

0.0015 

0.0016 

46 

0.3597 

0.4225 

1.1746 

0.0015 

0.0016 

48 

0.3622 

0.4276 

1.1806 

0.0015 

0.0016 

50 

0.3670 

0.4308 

1.1739 

0.0015 

0.0016 

Table  B.14:  Simulation  Results  for  Uniform  Populations  with  6  =  1.2  and  k  =  5 
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V 

PCSbem 

PCSavc 

Ratio 

s.e.  PCSbem 

s.e.  PCSavc 

2 

0.5003 

0.5866 

1.1725 

0.0016 

0.0016 

4 

0.5936 

0.6521 

1.0986 

0.0016 

0.0015 

6 

0.6460 

0.7140 

1.1052 

0.0015 

0.0014 

8 

0.6870 

0.7626 

1.1100 

0.0015 

0.0013 

10 

0.7266 

0.8000 

1.1009 

0.0014 

0.0013 

12 

0.7575 

0.8316 

1.0978 

0.0014 

0.0012 

14 

0.7826 

0.8562 

1.0940 

0.0013 

0.0011 

16 

0.8079 

0.8787 

1.0876 

0.0012 

0.0010 

18 

0.8274 

0.8960 

1.0830 

0.0012 

0.0010 

20 

0.8447 

0.9104 

1.0778 

0.0011 

0.0009 

22 

0.8618 

0.9245 

1.0728 

0.0011 

0.0008 

24 

0.8749 

0.9340 

1.0675 

0.0010 

0.0008 

26 

0.8874 

0.9436 

1.0633 

0.0010 

0.0007 

28 

0.8995 

0.9515 

1.0579 

0.0010 

0.0007 

30 

0.9080 

0.9576 

1.0546 

0.0009 

0.0006 

32 

0.9174 

0.9633 

1.0500 

0.0009 

0.0006 

34 

0.9256 

0.9680 

1.0459 

0.0008 

0.0006 

36 

0.9331 

0.9725 

1.0423 

0.0008 

0.0005 

38 

0.9388 

0.9760 

1.0396 

0.0008 

0.0005 

40 

0.9445 

0.9795 

1.0370 

0.0007 

0.0004 

42 

0.9499 

0.9819 

1.0336 

0.0007 

0.0004 

44 

0.9542 

0.9846 

1.0319 

0.0007 

0.0004 

46 

0.9582 

0.9867 

1.0298 

0.0006 

0.0004 

48 

0.9620 

0.9880 

1.0271 

0.0006 

0.0003 

50 

0.9660 

0.9894 

1.0243 

0.0006 

0.0003 

Table  B. 16:  Simulation  Results  for  Uniform  Populations  with  6  =  2.0  and  k 
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V 

PCSbem 

PCSavc 

Ratio 

s.e.  PCSbem 

s.e.  PCSavc 

2 

0.4009 

0.4967 

1.2391 

0.0016 

0.0016 

4 

0.4875 

0.5469 

1.1218 

0.0016 

0.0016 

6 

0.5323 

0.6157 

1.1566 

0.0016 

0.0015 

8 

0.5788 

0.6658 

1.1503 

0.0016 

0.0015 

0.6154 

0.7075 

1.1496 

0.0015 

0.0014 

12 

0.6515 

0.7446 

1.1429 

0.0015 

0.0014 

14 

0.6799 

0.7760 

1.1414 

0.0015 

0.0013 

16 

0.7080 

0.8028 

1.1339 

0.0014 

0.0013 

18 

0.7314 

0.8259 

1.1291 

0.0014 

0.0012 

0.7532 

0.8466 

1.1241 

0.0014 

0.0011 

22 

0.7755 

0.8650 

1.1154 

0.0013 

0.0011 

24 

0.7912 

0.8790 

1.1110 

0.0013 

0.0010 

26 

0.8083 

0.8920 

1.1035 

0.0012 

0.0010 

28 

0.8217 

0.9054 

1.1019 

0.0012 

0.0009 

0.8359 

0.9141 

1.0935 

0.0012 

0.0009 

32 

0.8483 

0.9242 

1.0894 

0.0011 

0.0008 

34 

0.8604 

0.9324 

1.0836 

0.0011 

0.0008 

36 

0.8708 

0.9396 

1.0790 

0.0011 

0.0008 

38 

0.8800 

0.9457 

1.0747 

0.0010 

0.0007 

0.8892 

0.9526 

1.0713 

0.0010 

0.0007 

42 

0.8984 

0.9570 

1.0652 

0.0010 

0.0006 

44 

0.9060 

0.9620 

1.0619 

0.0009 

0.0006 

46 

0.9123 

0.9659 

1.0588 

0.0009 

0.0006 

48 

0.9184 

0.9683 

1.0543 

0.0009 

0.0006 

0.9243 

0.9731 

1.0527 

0.0008 

0.0005 

Table  B.17:  Simulation  Results  for  Uniform  Populations  with  9  =  2.0  and  k  =  4 
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V 

PCSbem 

PCSavc 

Ratio 

s.e.  PCSbem 

s.e.  PCSavc 

2 

0.3337 

0.4289 

1.2852 

0.0015 

0.0016 

4 

0.4081 

0.4722 

1.1570 

0.0016 

0.0016 

6 

0.4513 

0.5375 

1.1911 

0.0016 

0.0016 

8 

0.4940 

0.5882 

1.1907 

0.0016 

0.0016 

10 

0.5293 

0.6304 

1.1910 

0.0016 

0.0015 

12 

0.5634 

0.6688 

1.1871 

0.0016 

0.0015 

14 

0.5924 

0.7028 

1.1863 

0.0016 

0.0014 

16 

0.6207 

0.7317 

1.1789 

0.0015 

0.0014 

18 

0.6457 

0.7576 

1.1732 

0.0015 

0.0014 

20 

0.6689 

0.7801 

1.1662 

0.0015 

0.0013 

22 

0.6921 

0.8025 

1.1596 

0.0015 

0.0013 

24 

0.7105 

0.8196 

1.1536 

0.0014 

0.0012 

26 

0.7294 

0.8353 

1.1452 

0.0014 

0.0012 

28 

0.7467 

0.8519 

1.1409 

0.0014 

0.0011 

30 

0.7610 

0.8643 

1.1358 

0.0013 

0.0011 

32 

0.7766 

0.8774 

1.1298 

0.0013 

0.0010 

34 

0.7905 

0.8872 

1.1223 

0.0013 

0.0010 

36 

0.8043 

0.8979 

1.1165 

0.0013 

0.0010 

38 

0.8158 

0.9073 

1.1122 

0.0012 

0.0009 

40 

0.8267 

0.9150 

1.1069 

0.0012 

0.0009 

42 

0.8368 

0.9235 

1.1036 

0.0012 

0.0008 

44 

0.8476 

0.9305 

1.0978 

0.0011 

0.0008 

46 

0.8564 

0.9365 

1.0935 

0.0011 

0.0008 

48 

0.8650 

0.9412 

1.0880 

0.0011 

0.0007 

50 

0.8732 

0.9471 

1.0847 

0.0011 

0.0007 

Table  B.18:  Simulation  Results  for  Uniform  Populations  with  0  =  2.0  and  k  =  5 


V 

PCSbem 

PCSavc 

Ratio 

s.e.  PCSbem 

s.e.  PCSavc 

2 

0.5438 

0.5533 

1.0175 

0.0016 

0.0016 

4 

0.5660 

0.5762 

1.0179 

0.0016 

0.0016 

6 

0.5826 

0.5940 

1.0196 

0.0016 

0.0016 

8 

0.5957 

0.6098 

1.0237 

0.0016 

0.0015 

m 

0.6082 

0.6238 

1.0257 

0.0015 

0.0015 

12 

0.6181 

0.6351 

1.0275 

0.0015 

0.0015 

14 

0.6285 

0.6461 

1.0281 

0.0015 

0.0015 

16 

0.6379 

0.6565 

1.0292 

0.0015 

0.0015 

18 

0.6458 

0.6662 

1.0316 

0.0015 

0.0015 

0.6536 

0.6732 

1.0300 

0.0015 

0.0015 

22 

0.6615 

0.6813 

1.0299 

0.0015 

0.0015 

24 

0.6678 

0.6888 

1.0315 

0.0015 

0.0015 

26 

0.6745 

0.6957 

1.0314 

0.0015 

0.0015 

28 

0.6802 

0.7030 

1.0335 

0.0015 

0.0014 

0.6862 

0.7088 

1.0329 

0.0015 

0.0014 

32 

0.6923 

0.7145 

1.0321 

0.0015 

0.0014 

34 

0.6980 

0.7201 

1.0316 

0.0015 

0.0014 

36 

0.7039 

0.7258 

1.0312 

0.0014 

0.0014 

38 

0.7092 

0.7315 

1.0313 

0.0014 

0.0014 

0.7146 

0.7374 

1.0319 

0.0014 

0.0014 

42 

0.7201 

0.7430 

1.0318 

0.0014 

0.0014 

44 

0.7227 

0.7481 

1.0351 

0.0014 

0.0014 

46 

0.7276 

0.7530 

1.0349 

0.0014 

0.0014 

48 

0.7325 

0.7572 

1.0337 

0.0014 

0.0014 

0.7371 

0.7613 

1.0328 

0.0014 

0.0013 

Table  B.19:  Simulation  Results  for  Gamma  Populations  with  6  =  1.2  and  k 


V 

PCSbem 

PCSavc 

Ratio 

s.e.  PCSbem 

s.e.  PCSavc 

2 

0.3742 

0.3863 

1.0323 

0.0015 

0.0015 

4 

0.3974 

0.4071 

1.0243 

0.0015 

0.0016 

6 

0.4112 

0.4241 

1.0313 

0.0016 

0.0016 

8 

0.4211 

0.4392 

1.0429 

0.0016 

0.0016 

10 

0.4340 

0.4531 

1.0441 

0.0016 

0.0016 

12 

0.4420 

0.4648 

1.0515 

0.0016 

0.0016 

14 

0.4506 

0.4752 

1.0547 

0.0016 

0.0016 

16 

0.4608 

0.4853 

1.0532 

0.0016 

0.0016 

18 

0.4689 

0.4962 

1.0584 

0.0016 

0.0016 

20 

0.4740 

0.5048 

1.0650 

0.0016 

0.0016 

22 

0.4832 

0.5132 

1.0621 

0.0016 

0.0016 

24 

0.4898 

0.5222 

1.0662 

0.0016 

0.0016 

26 

0.4958 

0.5296 

1.0682 

0.0016 

0.0016 

28 

0.5020 

0.5369 

1.0696 

0.0016 

0.0016 

30 

0.5076 

0.5432 

1.0700 

0.0016 

0.0016 

32 

0.5140 

0.5500 

1.0700 

0.0016 

0.0016 

34 

0.5199 

0.5567 

1.0708 

0.0016 

0.0016 

36 

0.5256 

0.5627 

1.0707 

0.0016 

0.0016 

38 

0.5309 

0.5693 

1.0723 

0.0016 

0.0016 

40 

0.5375 

0.5758 

1.0712 

0.0016 

0.0016 

42 

0.5435 

0.5828 

1.0722 

0.0016 

0.0016 

44 

0.5481 

0.5882 

1.0732 

0.0016 

0.0016 

46 

0.5538 

0.5944 

1.0732 

0.0016 

0.0016 

48 

0.5585 

0.5995 

1.0734 

0.0016 

0.0016 

50 

0.5628 

0.6047 

1.0744 

0.0016 

0.0015 

Table  B.20:  Simulation  Results  for  Gamma  Populations  with  9  =  1.2  and  k 
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V 

PCSbem 

PCSavc 

Ratio 

s.e.  PCSbem 

s.e.  PCSavc 

2 

0.2863 

0.2962 

1.0346 

0.0014 

0.0014 

4 

0.3069 

0.3180 

1.0361 

0.0015 

0.0015 

6 

0.3153 

0.3328 

1.0558 

0.0015 

0.0015 

8 

0.3263 

0.3445 

1.0557 

0.0015 

0.0015 

10 

0.3355 

0.3559 

1.0607 

0.0015 

0.0015 

12 

0.3431 

0.3687 

1.0745 

0.0015 

0.0015 

14 

0.3505 

0.3758 

1.0722 

0.0015 

0.0015 

16 

0.3573 

0.3847 

1.0765 

0.0015 

0.0015 

18 

0.3637 

0.3930 

1.0806 

0.0015 

0.0015 

20 

0.3696 

0.4008 

1.0844 

0.0015 

0.0016 

22 

0.3755 

0.4068 

1.0832 

0.0015 

0.0016 

24 

0.3805 

0.4145 

1.0896 

0.0015 

0.0016 

26 

0.3862 

0.4218 

1.0922 

0.0015 

0.0016 

28 

0.3924 

0.4283 

1.0916 

0.0015 

0.0016 

30 

0.3977 

0.4354 

1.0947 

0.0015 

0.0016 

32 

0.4028 

0.4419 

1.0970 

0.0016 

0.0016 

34 

0.4081 

0.4481 

1.0980 

0.0016 

0.0016 

36 

0.4129 

0.4538 

1.0992 

0.0016 

0.0016 

38 

0.4170 

0.4598 

1.1027 

0.0016 

0.0016 

40 

0.4225 

0.4654 

1.1016 

0.0016 

0.0016 

42 

0.4265 

0.4710 

1.1045 

0.0016 

0.0016 

44 

0.4312 

0.4773 

1.1071 

0.0016 

0.0016 

46 

0.4357 

0.4831 

1.1089 

0.0016 

0.0016 

48 

0.4406 

0.4888 

1.1093 

0.0016 

0.0016 

50 

0.4445 

0.4946 

1.1126 

0.0016 

0.0016 

Table  B.21:  Simulation  Results  for  Gamma  Populations  with  6  =  1.2  and  k  =  4 
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V 

PCSbem 

PCSavc 

Ratio 

s.e.  PCSbem 

s.e.  PCSavc 

2 

0.2393 

1.0337 

0.0013 

0.0013 

4 

0.2589 

1.0419 

0.0014 

0.0014 

6 

0.2705 

1.0565 

0.0014 

0.0014 

8 

0.2801 

1.0579 

0.0014 

0.0014 

10 

0.2907 

1.0630 

0.0014 

0.0014 

12 

0.2802 

0.3025 

1.0794 

0.0014 

0.0015 

14 

0.2850 

0.3094 

1.0857 

0.0014 

0.0015 

16 

0.2914 

0.3166 

1.0862 

0.0014 

0.0015 

18 

0.2973 

0.3246 

1.0919 

0.0014 

0.0015 

20 

0.3022 

0.3308 

1.0947 

0.0015 

0.0015 

0.3365 

1.0953 

0.0015 

0.0015 

0.3431 

1.1014 

0.0015 

0.0015 

26 

0.3168 

0.3492 

1.1024 

0.0015 

0.0015 

28 

0.3214 

0.3548 

1.1041 

0.0015 

0.0015 

30 

0.3268 

0.3612 

1.1053 

0.0015 

0.0015 

32 

0.3298 

0.3667 

1.1119 

0.0015 

0.0015 

34 

0.3337 

0.3732 

1.1183 

0.0015 

0.0015 

36 

0.3380 

0.3796 

1.1231 

0.0015 

0.0015 

38 

0.3424 

0.3844 

1.1227 

0.0015 

0.0015 

40 

0.3464 

0.3893 

1.1240 

0.0015 

0.0015 

42 

0.3516 

0.3948 

1.1228 

0.0015 

0.0015 

44 

0.3551 

0.4001 

1.1267 

0.0015 

0.0015 

46 

0.3586 

0.4054 

1.1307 

0.0015 

0.0016 

48 

0.3625 

0.4103 

1.1318 

0.0015 

0.0016 

50 

0.3665 

0.4154 

1.1334 

0.0015 

0.0016 

Table  B.22:  Simulation  Results  for  Gamma  Populations  with  6  =  1.2  and  k 
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V 

pCSbem 

PCSavc 

Ratio 

s.e.  PCSbem 

s.e.  PCSavc 

2 

0.6672 

0.7005 

HI 

0.0015 

0.0014 

4 

0.7405 

0.7782 

EEH 

0.0014 

0.0013 

6 

0.7895 

0.8267 

1.0472 

0.0013 

0.0012 

8 

0.8254 

0.8608 

1.0429 

0.0012 

0.0011 

10 

0.8553 

0.8871 

1.0372 

0.0011 

0.0010 

12 

0.8782 

0.9077 

1.0336 

0.0010 

0.0009 

14 

0.8959 

0.9242 

1.0315 

0.0010 

0.0008 

16 

0.9120 

0.9377 

1.0282 

0.0009 

0.0008 

18 

0.9248 

0.9481 

1.0252 

0.0008 

0.0007 

20 

0.9353 

0.9569 

1.0230 

0.0008 

0.0006 

22 

0.9446 

0.9636 

1.0202 

0.0007 

0.0006 

24 

0.9517 

0.9700 

1.0192 

0.0007 

0.0005 

26 

0.9578 

0.9747 

1.0177 

0.0006 

0.0005 

28 

0.9633 

0.9785 

1.0157 

0.0006 

0.0005 

30 

0.9681 

0.9819 

1.0142 

0.0006 

0.0004 

32 

0.9722 

0.9846 

1.0128 

0.0005 

0.0004 

34 

0.9760 

0.9872 

1.0115 

0.0005 

0.0004 

36 

0.9790 

0.9890 

1.0102 

0.0005 

0.0003 

38 

0.9815 

0.9905 

1.0091 

0.0004 

0.0003 

40 

0.9842 

0.9917 

1.0077 

0.0004 

0.0003 

42 

0.9861 

0.9930 

1.0070 

0.0004 

0.0003 

44 

0.9878 

0.9940 

1.0063 

0.0003 

0.0002 

46 

0.9893 

0.9948 

1.0056 

0.0003 

0.0002 

48 

0.9955 

1.0050 

0.0003 

0.0002 

0.9961 

1.0044 

0.0003 

0.0002 

Table  B.23:  Simulation  Results  for  Gamma  Populations  with  9  =  2.0  and  k 


V 

PCSbem 

PCSavc 

Ratio 

s.e.  PCSbem 

s.e.  PCSavc 

2 

0.4008 

0.4479 

1.1174 

0.0016 

0.0016 

4 

0.4874 

0.5372 

1.1023 

0.0016 

0.0016 

6 

0.5324 

0.5993 

1.1256 

0.0016 

0.0016 

8 

0.5803 

0.6506 

1.1211 

0.0016 

0.0015 

0.6175 

0.6945 

1.1247 

0.0015 

0.0015 

12 

0.6508 

0.7311 

1.1233 

0.0015 

0.0014 

14 

0.6800 

0.7624 

1.1211 

0.0015 

0.0013 

16 

0.7901 

1.1172 

0.0014 

0.0013 

18 

0.7317 

0.8145 

1.1132 

0.0014 

0.0012 

0.7528 

0.8356 

1.1099 

0.0014 

0.0012 

22 

0.7729 

0.8536 

1.1045 

0.0013 

0.0011 

24 

Essa 

0.8702 

1.1013 

0.0013 

0.0011 

26 

m 

0.8849 

1.0973 

0.0012 

0.0010 

28 

0.8214 

0.8973 

1.0924 

0.0012 

0.0010 

0.8354 

0.9085 

1.0876 

0.0012 

0.0009 

32 

0.9183 

1.0842 

0.0011 

0.0009 

34 

0.8589 

0.9265 

1.0787 

0.0011 

0.0008 

36 

0.8696 

0.9349 

1.0752 

0.0011 

0.0008 

38 

0.8786 

0.9417 

1.0718 

0.0010 

0.0007 

0.8891 

0.9476 

1.0658 

0.0010 

0.0007 

42 

0.8971 

0.9533 

1.0626 

0.0010 

0.0007 

44 

0.9582 

1.0587 

0.0009 

0.0006 

46 

0.9121 

0.9626 

1.0554 

0.0009 

0.0006 

48 

0.9660 

1.0513 

0.0009 

0.0006 

m 

0.9246 

0.9696 

1.0486 

0.0008 

0.0005 

Table  B.25:  Simulation  Results  for  Gamma  Populations  with  6  =  2.0  and  k  =  4 
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V 

PCSbem 

PCSavc 

Ratio 

s.e.  PCSbem 

s.e.  PCSavc 

2 

0.3334 

0.3758 

1.1273 

0.0015 

0.0015 

4 

0.4086 

0.4594 

1.1243 

0.0016 

0.0016 

6 

0.4519 

0.5211 

1.1530 

0.0016 

0.0016 

8 

0.4946 

0.5728 

1.1581 

0.0016 

0.0016 

10 

0.5330 

0.6178 

1.1591 

0.0016 

0.0015 

12 

0.5659 

0.6562 

1.1595 

0.0016 

0.0015 

14 

0.5945 

0.6898 

1.1605 

0.0016 

0.0015 

16 

0.6210 

0.7189 

1.1578 

0.0015 

0.0014 

18 

0.6477 

0.7469 

1.1533 

0.0015 

0.0014 

20 

0.6702 

0.7708 

1.1501 

0.0015 

0.0013 

22 

0.6913 

0.7915 

1.1449 

0.0015 

0.0013 

24 

0.7105 

0.8124 

1.1434 

0.0014 

0.0012 

26 

0.7294 

0.8297 

1.1375 

0.0014 

0.0012 

28 

0.7457 

0.8460 

1.1344 

0.0014 

0.0011 

30 

0.7617 

0.8606 

1.1298 

0.0013 

0.0011 

32 

0.7757 

0.8726 

1.1249 

0.0013 

0.0011 

34 

0.7900 

0.8841 

1.1190 

0.0013 

0.0010 

36 

0.8037 

0.8946 

1.1131 

0.0013 

0.0010 

38 

0.8151 

0.9041 

1.1092 

0.0012 

0.0009 

40 

0.8259 

0.9130 

1.1054 

0.0012 

0.0009 

42 

0.8370 

0.9204 

1.0995 

0.0012 

0.0009 

44 

0.8474 

0.9275 

1.0946 

0.0011 

0.0008 

46 

0.8567 

0.9340 

1.0902 

0.0011 

0.0008 

48 

0.8658 

0.9399 

1.0857 

0.0011 

0.0008 

50 

0.8733 

0.9449 

1.0820 

0.0011 

0.0007 

Table  B.26:  Simulation  Results  for  Gamma  Populations  with  9  =  2.0  and  k 


V 

PCSbem 

PCSavc 

Ratio 

s.e.  PCSbem 

s.e.  PCSavc 

2 

0.5453 

0.5566 

1.0208 

0.0016 

0.0016 

4 

0.5671 

0.5920 

1.0438 

0.0016 

0.0016 

6 

0.5838 

0.6148 

1.0532 

0.0016 

0.0015 

8 

0.5988 

0.6357 

1.0615 

0.0016 

0.0015 

10 

0.6100 

0.6514 

1.0679 

0.0015 

0.0015 

12 

0.6215 

0.6661 

1.0717 

0.0015 

0.0015 

14 

0.6309 

0.6782 

1.0749 

0.0015 

0.0015 

16 

0.6395 

0.6901 

1.0790 

0.0015 

0.0015 

18 

0.6479 

0.7034 

1.0857 

0.0015 

0.0014 

20 

0.6564 

0.7132 

1.0866 

0.0015 

0.0014 

22 

0.6631 

0.7224 

1.0894 

0.0015 

0.0014 

24 

0.6697 

0.7320 

1.0932 

0.0015 

0.0014 

26 

0.6764 

0.7408 

1.0952 

0.0015 

0.0014 

28 

0.6839 

0.7503 

1.0971 

0.0015 

0.0014 

30 

0.6896 

0.7576 

1.0986 

0.0015 

0.0014 

32 

0.6960 

0.7642 

1.0981 

0.0015 

0.0013 

34 

0.7018 

0.7716 

1.0996 

0.0014 

0.0013 

36 

0.7070 

0.7785 

1.1013 

0.0014 

0.0013 

38 

0.7118 

0.7854 

1.1034 

0.0014 

0.0013 

40 

0.7159 

0.7907 

1.1046 

0.0014 

0.0013 

42 

0.7209 

0.7955 

1.1036 

0.0014 

0.0013 

44 

0.7250 

0.8016 

1.1056 

0.0014 

0.0013 

46 

0.7293 

0.8073 

1.1070 

0.0014 

0.0012 

48 

0.7332 

0.8122 

1.1076 

0.0014 

0.0012 

50 

0.7387 

0.8179 

1.1072 

0.0014 

0.0012 

Table  B.27:  Simulation  Results  for  Bernoulli  Populations  with  6  =  1.2  and  k 
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V 

PCSbem 

PCSavc 

Ratio 

s.e.  PCSbem 

s.e.  PCSavc 

2 

0.3744 

0.3908 

1.0438 

0.0015 

0.0015 

4 

0.3996 

0.4244 

1.0621 

0.0015 

0.0016 

6 

0.4126 

0.4435 

1.0748 

0.0016 

0.0016 

8 

0.4240 

0.4647 

1.0961 

0.0016 

0.0016 

10 

0.4357 

0.4804 

1.1026 

0.0016 

0.0016 

12 

0.4446 

0.4948 

1.1128 

0.0016 

0.0016 

14 

0.4523 

0.5069 

1.1208 

0.0016 

0.0016 

16 

0.4616 

0.5204 

1.1275 

0.0016 

0.0016 

18 

0.4690 

0.5317 

1.1337 

0.0016 

0.0016 

20 

0.4770 

0.5442 

1.1407 

0.0016 

0.0016 

22 

0.4841 

0.5525 

1.1413 

0.0016 

0.0016 

24 

0.4910 

0.5656 

1.1519 

0.0016 

0.0016 

26 

0.4984 

0.5740 

1.1518 

0.0016 

0.0016 

28 

0.5041 

0.5845 

1.1595 

0.0016 

0.0016 

30 

0.5106 

0.5930 

1.1613 

0.0016 

0.0016 

0.5161 

0.6017 

1.1658 

0.0016 

0.0015 

0.6117 

1.1711 

0.0016 

0.0015 

36 

0.5279 

0.6196 

1.1738 

0.0016 

0.0015 

38 

0.5338 

0.6277 

1.1759 

0.0016 

0.0015 

40 

0.5382 

0.6337 

1.1774 

0.0016 

0.0015 

42 

0.5436 

0.6408 

1.1788 

0.0016 

0.0015 

44 

0.5487 

0.6483 

1.1814 

0.0016 

0.0015 

46 

0.5536 

0.6564 

1.1857 

0.0016 

0.0015 

48 

0.5581 

0.6623 

1.1868 

0.0016 

0.0015 

50 

0.5635 

0.6698 

1.1885 

0.0016 

0.0015 

Table  B.28:  Simulation  Results  for  Bernoulli  Populations  with  6  =  1.2  and  k  =  3 
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Appendix  C 

PROGRAM  LISTINGS 


C**************************************************************** 
c************************************* ******* ******************** 
£**************************************************************** 
c 

c  THIS  IS  PROGRAM  FOR  MULTINOMIAL  SELECTION  PROBLEM.  METHOD  I 
c  IS  BECHHOFER’S  METHOD  (BEM) .  METHOD  II  IS  OUR  NEW  APPROACH 
c  MAKING  ALL  VECTOR  COMPARISONS  (AVC) .  BEST  POPULATION  IS 
c  MODELED  AS  XI.  SETUP  TO  USE  THE  FOLLOWING  DISTRIBUTIONS: 
c 

C  EXPONENTIAL  -  SET  LAMBDA=1.0,  VARY  MU  TO  MAINTAIN  THETA 

C  CONTINUOUS  UNIFORM  -  SET  A=l,  VARY  B  TO  MAINTAIN  THETA 

C  GAMMA  -  SET  A1=A2=3.0,  B2=1.0,  VARY  B1  TO  MAINTAIN  THETA 

c 

c  THIS  IS  REVISED  VERSION  SORTING  AND  THEN  COMPARING  EACH  VALUE 
c  AGAINST  ALL  VALUES  FROM  REMAINING  POPULATIONS  INSTEAD  OF 
c  ACTUALLY  FORMING  ALL  V~K  PSEUDO-REPLICATIONS, 
c 

c  SETUP  FOR  MAX  OF  5  POPULATIONS  AND  100  VECTORS, 
c 

£***************************************************************** 
£***************************************************************** 
c  VARIABLE  TABLE 

c 

c  NAME  DESCRIPTION 

c 

c  NREPS  NUMBER  OF  REPLICATIONS 

c 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


NVECT  NUMBER  OF  INDEPENDENT  VECTORS  TO  SIMULATE. 

SET  UP  FOR  MAXIMUM  OF  100. 

NPOP  NUMBER  OF  POPULATIONS.  SET  UP  FOR  MAXIMUM 

OF  5. 

LAMBDA  EXPONENTIAL  RATE  FOR  POPULATION  1,  WHICH  IS 

ALWAYS  MODELED  AS  BEST  SYSTEM.  SET  TO  1.0. 

MU  EXPONENTIAL  RATE  FOR  POPULATIONS  OTHER  THAN 

1.  USE  OF  A  SINGLE  MU  ASSUMES  METHOD  I 
l.f.c.  VARIED  WITH  INCREASING  K  TO  MAINTAIN 
DESIRED  THETA  =  P(X>Y)/(1-P(X>Y)/K-1) . 

A  RIGHT  END  POINT  FOR  CONTINUOUS  UNIFORM  DIST 

FOR  POPULATION  1.  SET  TO  1.0. 

B  RIGHT  END  POINT  FOR  CONTINUOUS  UNIFORM  DIST 

FOR  POPULATIONS  OTHER  THAN  1.  VARIED  WITH 
INCREASING  K  TO  MAINTAIN  DESIRED  THETA. 


JX1(3)-JX5(3) 

RUNF(JX) 

PCS1(V) 

PCS2(V) 

SE1(V) 

SE2(V) 

RTO(V) 

NWIN1 (L ,M) 

X( J,K) 


RANDOM  NUMBER  SEEDS  FOR  USE  WITH  L’ECUYER’S 
GENERATOR.  SEPARATE  STREAM  FOR  EACH  POP. 

CALL  TO  L’ECUYER’S  GENERATOR. 

PCS  FOR  METHOD  I  FOR  V  =  THE  NUMBER  OF 
VECTORS,  V=l,2, . . . ,100. 

PCS  FOR  METHOD  II  FOR  V. 

STANDARD  ERROR  FOR  PCS1(V). 

STANDARD  ERROR  FOR  PCS2(V). 

RATIO  OF  PCS2(V)/PCS1 (V) . 

NUMBER  OF  WINS  WITH  METHOD  I  FOR  POPULATION 
L  WITH  M  VECTORS;  L=2,3,...,K;  M=l,2, . . . ,100. 

OBSERVATION  K  FROM  POPULATION  J. 
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USED  TO  KEEP  TRACK  OF  HOW  MANY  OBSERVATIONS 
FROM  ALL  OTHER  POPULATIONS  AN  OBSERVATION 
FROM  POPULATION  J  EXCEEDS. 

USED  TO  STORE  HOW  MANY  OBSERVATIONS  IN  POP  M 
ARE  EXCEEDED  BY  CURRENT  OBSERVATION  FROM 
POP  J.  PRODUCT  OF  XJ(J,M)  FOR  M=2,...,K; 

M  NOT  EQUAL  TO  K  RESULTS  IN  TOTAL  XCNT(J) 

FOR  A  SINGLE  OBSERVATION  FROM  POP  J. 

SUM  OF  XCNT(J)  OVER  ALL  OBSERVATIONS  FOR 
POP  J. 
c 

c 

PROGRAM  MSP 
c 
c 

REAL  LAMBDA, MU, X(5, 100) ,XJ(5,5) ,XCNT(5) ,FC(5) 

REAL  PCSl(lOO) ,PCS2(100) ,SE1(100) ,SE2(100) ,RT0(100) 

INTEGER* 4  JX1(3) , JX2(3) , JX3(3) , JX4(3) , JX5(3) 

INTEGER* 4  JX6(3) , JX7(3) , JX8(3) , JX9(3) , JX10(3) 

INTEGER  NWIN1 (5 , 100) 
c 

c  FILES  TO  OUTPUT  RESULTS  TO 

c 

OPEN (4 , FILE= ’ exk5t20 . dat ’ , STATUS= ’ UNKNOWN » ) 
c 

c  INITIALIZE  VARIABLES 

c 

c  IDIST=1  IF  EXPONENTIAL,  2  IF  CONTINUOUS  UNIFORM 

c  3  IF  GAMMA 

c 

NREPS=100000 
NVECT=2 
NP0P=5 
LAMBDA=1 . 0 
MU= 1.5885 
A=1 . 0 
B=1 . 1 
AX=3 . 0 


c  XCNT(J) 

c 
c 
c 

c  XJ(J,M) 

c 

c 

c 

c 

c 

c  FC( J) 

c 
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c 


BX=1 .3061 
AY=3 . 0 
BY=1 . 0 
IDIST=1 

c  SETUP  CONSTANTS  FOR  GAMMA  RV  GENERATION 
c 

AB=1 ,0/SQRT(2*AX-l) 

AR=1.0/SQRT(2*AY-1) 

BB=AX-L0G(4 . 0) 

BR=AY-L0G(4.0) 

QB=AX+1 .O/AB 
QR=AY+1 . O/AR 
TB=4 . 5 
TR=4 . 5 

DB=1 . O+L0G(TB) 

DR=1 . O+LOG (TR) 
c 

c  SEEDS  FOR  L’ECUYER’S  RN  GENERATOR 
c 

JX1 (2)=748932582 
JX1 (3)=639287 
JX2 (2) =64298628 
JX2(3)=196998 
JX3 (2) =40689408 
JX3 (3) =122595154 
JX4 (2) =680620100 
JX4 (3) =3350831 18 
JX5 (2) =64918046 
JX5 (3) =812715188 
c 

5  DO  10  1=1 ,NVECT 
PCS1(I)=0 
PCS2(I)=0 
10  CONTINUE 

DO  2000  1=1 ,NREPS 
DO  30  J=1 ,NPOP 
FC( J)=0 

DO  20  K=1 ,NVECT 
NWIN1 (J,K)=0 

20  CONTINUE 
30  CONTINUE 
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c  RANDOM  VARIATE  GENERATION  DONE  HERE.  SEPARATE  SET  OF  RN 
c  STREAMS  FOR  EACH  POPULTION. 
c 

IF(IDIST.EQ.2)G0  TO  65 
IF(IDIST.EQ.3)G0  TO  81 
c 

c  EXPONENTIAL  RANDOM  VARIATES 
c 

DO  60  K=1 ,NVECT 

X ( 1 , K) =- (1/LAMBDA) *LOG ( 1-RUNF ( JX1 ) ) 

X(2 ,K)=-(1/MU) *L0G(1-RUNF( JX2) ) 

IF (NPOP .EQ . 2)G0  TO  60 
X (3 , K) =-( 1/MU) *LOG ( 1-RUNF ( JX3) ) 

IF (NPOP. EQ. 3) GO  TO  60 
X (4 , K) =-( 1/MU) *LOG ( 1-RUNF ( JX4) ) 

IF (NPOP. EQ. 4) GO  TO  60 

X(5 ,K)=-(1/MU) *LOG ( 1-RUNF (JX5) ) 

60  CONTINUE 

GO  TO  89 
c 

c  CONTINOUS  UNIFORM  RANDOM  VARIATES 
c 

65  DO  80  K=1 ,NVECT 

X(1 ,K)=B*RUNF ( JX1) 

X(2,K)=A*RUNF(JX2) 

IF (NPOP. EQ. 2) GO  TO  80 
X(3,K)=A*RUNF( JX3) 

IF (NPOP. EQ. 3) GO  TO  80 
X(4 ,K)=A*RUNF( JX4) 

IF (NPOP. EQ. 4) GO  TO  80 
X(5 ,K)=A*RUNF( JX5) 

80  CONTINUE 

GO  TO  89 
c 

c  GAMMA  RANDOM  VARIATES 
c 

81  DO  88  K=1 ,NVECT 

82  U1=RUNF( JX1) 

U2=RUNF ( JX2) 

V=AB*L0G(U1/ (1-U1) ) 

Y=AX*EXP(V) 
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Z=(U1**2)*U2 

W=BB+QB*V-Y 

IF (W+DB-TB*Z . GE . 0) THEN 
XT=Y 

ELSE  IF (W . GE . LOG (Z) ) THEN 
XT=Y 
ELSE 

GO  TO  82 
END  IF 

X(1 ,K)=BX*XT 

83  U1=RUNF( JX3) 

U2=RUNF ( JX4) 

V=AR*L0G(U1/ (1-U1)) 
Y=AY*EXP(V) 

Z=(U1**2)*U2 

W=BR+QR*V-Y 

IF (W+DR-TR*Z . GE . 0) THEN 
OT=Y 

ELSE  IF (W . GE . LOG (Z) ) THEN 
OT=Y 
ELSE 

GO  TO  83 
END  IF 

X(2,K)=BY*0T 
IF(NP0P.EQ.2)G0  TO  88 

84  U1=RUNF( JX5) 

U2=RUNF( JX6) 
V=AR*L0G(U1/(1-U1) ) 
Y=AY*EXP(V) 

Z= (Ul**2) *U2 
W=BR+QR*V-Y 
IF (W+DR-TR*Z . GE . 0) THEN 
OT=Y 

ELSE  IF (W . GE . LOG (Z) ) THEN 
OT=Y 
ELSE 

GO  TO  84 
END  IF 

X(3,K)=BY*0T 
IF(NP0P.EQ.3)G0  TO  88 

85  U1=RUNF( JX7) 

U2=RUNF(JX8) 
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V=AR*L0G(U1/(1-U1) ) 

Y=AY*EXP(V) 

Z=(U1**2)*U2 
W=BR+QR*V-Y 
I F ( W+DR-TR*  Z . GE . 0 ) THEN 
OT=Y 

ELSE  IF ( W . GE . LOG (Z) ) THEN 
0T=Y 
ELSE 

GO  TO  85 
END  IF 

X(4,K)=BY*OT 
IF(NP0P.EQ.4)G0  TO  88 
86  U1=RUNF ( JX9) 

U2=RUNF ( JXIO) 

V=AR*LQG (U1/C1-U1) ) 

Y=AY*EXP(V) 

Z=(U1**2) *U2 
W=BR+QR*V-Y 
IF ( W+DR-TR* Z . GE . 0 ) THEN 
OT=Y 

ELSE  IF (W . GE . LOG (Z) ) THEN 
OT=Y 
ELSE 

GO  TO  86 
END  IF 

X(5 ,K)=BY*OT 

88  CONTINUE 
c 

c  METHOD  I  (BEM) 

c  SIMPLY  NEED  TO  COUNT  NUMBER  OF  TIMES  EACH  Xi  IS  BEST  IN 

c  SAME  VECTOR, 

c 

89  DO  100  K=1 ,NVECT 

BEST=X(1 ,K) 

DO  90  J=2 ,NPOP 

TM= AM AX 1 ( BEST , X ( J , K ) ) 

IF (TM . GT . BEST) BEST=TM 

90  CONTINUE 

DO  95  L=1 ,NPOP 

IF(X(L,K) .EQ.TM)NWIN1(L,K)=NWIN1(L,K)+1 
95  CONTINUE 
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100 


CONTINUE 


c 

c  IF  BEST  POPULATION  (#1)  WINS  OUTRIGHT,  INCREASE  PCS  BY  ONE. 
c  IF  ANOTHER  POPULATION  WINS,  DO  NOT  INCREASE  PCS. 
c  FOR  TIES  INCLUDING  BEST  POPULATION,  NEED  TO  KEEP  TRACK  OF  THE 

c  NUMBER  OF  TIES  AND  INCREASE  PCS  BY  1/NTIES. 

c 

DO  105  K=2 , NVECT 
DO  102  J=1 ,NPOP 

NWIN1 ( J ,K)=NWIN1 ( J,K)+NWIN1 ( J ,K-1) 

102  CONTINUE 
105  CONTINUE 

DO  111  K=l, NVECT 
NTIES=1 

DO  110  J=2 ,NPOP 

IF (NWIN1 ( 1 , K) . LT . NWIN1 ( J , K) ) THEN 
PCS1 (K)=PCS1 (K) 

GO  TO  111 
ELSE 

IF (NWIN1 ( 1 , K) . EQ . NWIN1 ( J ,K) ) NTIES=NTIES+1 
IF ( J . Eq . NPOP) PCS 1 (K) =PCS1 (K) +1 . 0/ (NTIES) 

END  IF 

110  CONTINUE 

111  CONTINUE 
c 

c  METHOD  II  (AVC) 

c 

DO  1180  NM2=1, NVECT 
CALL  SORTV (NPOP , NM2 , X) 

DO  1150  K=1 ,NM2 
DO  1148  J=1 ,NPOP 
XCNT(J)=1 
DO  1146  M=1 ,NPOP 
DO  1142  L=1 ,NM2 

IFCM.NE. J. AND. X(J,K) . GT.X(M,L))THEN 
XJ( J,M)=NM2- (L-l) 

GO  TO  1145 
ELSE 

IF (L . EQ . NM2 )XJ(J,M)=0 
END  IF 

1142  CONTINUE 

1145  IFCM.NE. J)XCNT(J)=XCNT(J)*XJ(J,M) 
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1146  CONTINUE 

1148  CONTINUE 

DO  1149  J=1 ,NPOP 

FC( J)=FC( J)+XCNT( J) 

1149  CONTINUE 

1150  CONTINUE 
NT=1 

DO  1170  J=2 ,NPOP 

IF(FC(1) .LT.FC(J) )THEN 
PCS2 (NM2) =PCS2 (NM2) 

GO  TO  1180 
ELSE 

IF (FC ( 1 ) .EQ.FC(J))NT=NT+1 
IF  (  J .  EC) .  NPOP)  PCS2  (NM2)  =PCS2  (NM2)  +1 . 0/NT 
END  IF 

1170  CONTINUE 
1180  CONTINUE 
2000  CONTINUE 

DO  2010  1=1 ,NVECT 

PCS1(I)=PCS1(I) /REAL ( NREPS ) 

PCS2 ( I ) =PCS2 ( I ) /REAL (NREPS) 

SE1 (I)=(PCS1 (I)*(l . 0-PCS1 (I) ) /REAL (NREPS) )** . 5 
SE2(I)=(PCS2(I)*(1 .0-PCS2 (I)) /REAL (NREPS))** .5 
RT0(I)=PCS2(I)/PCS1 (I) 
c 

c  FILE  4  SETUP  TO  READ  DIRECTLY  INTO  MINITAB  OR  OTHER  PROGRAM 
c 

WRITE (4, 2050) NPOP, I, PCS 1(1) ,PCS2(I) ,RTO(I) ,SE1(I) ,SE2(I) 
c 

2010  CONTINUE 

2050  FORMAT ( IX , 214 , 2X , 5 (F8 . 6 , 2X) ) 

STOP 

END 

c 

SUBROUTINE  SORTV  (K,N, Y) 
c 

INTEGER  K,N 
REAL  Y (5 , 100) 
c 

DO  30  L=1,K 
DO  20  J=2,N 
TEMP=Y(L, J) 
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I=J-1 

10  IF(I.EQ.O)GO  TO  20 

IF (TEMP . GT . Y (L , I ) ) THEM 
Y(L,I+1)=Y(L,I) 

Y (L , I)=TEMP 
1=1-1 
GO  TO  10 
END  IF 
20  CONTINUE 
30  CONTINUE 
RETURN 
END 
c 

FUNCTION  RUNF(JX) 

c  L’Ecuyer’s  generator  as  given  as  unifl  in  BF&S,  p.  332 
c 

c  inputs: 

c  jx(2),  jx(3)  =  two  random  integers  with: 

c  0  <  jx(2)  <  2147483563 

c  0  <  jx(3)  <  2147483399 

c 

c  outputs: 

c  jx(l)  pseudorandom  integer,  0  <  jx(3)  <  2147483563 

c  period  of  jx(l)  is  about  2.3el8 

c  jx(2)  =  pseudorandom  integer,  0  <  jx(2)  <  2147483563 

c  jx(3)  =  pseudorandom  integer,  0  <  jx(2)  <  2147483399 

c  unifl  =  pseudorandom  real,  0.  <  unifl  <  1. 

c 

Q3|c3|ca|e9(c3|c3(c34c3fc9|e3|ea|c3fc9tca4c9|c3fc3fca|e^c3fe9|c3|e3fe3|cafc9|c34c34e3|e9fc34ca|c9|ea^c3|eaiea|c3|e3|ea|ca|c3|ca4c34c99ea|e3|c3fca|c9|c94c3fc9|e9|e3|c3|ca|c3|c9|ca4e9|c 

INTEGERS  JX,  K 
DIMENSION  JX(3) 
c 

c  get  next  term  in  the  first  stream  =  40014* jx(2)  mod  2147483563 
K  =  JX(2)/53668 

JX(2)  =  40014*(JX(2)  -  K*53668)  -  K*12211 
IF  ( JX(2)  .LT.  0)  JX(2)  =  JX(2)  +  2147483563 
c 

c  get  next  term  in  the  second  stream  =  40692*jx(3)  mod  2147483399 
K  =  JX(3)/52774 

JX(3)  =  40692* (JX(3)  -  K*52774)  -  K*3791 
IF  ( JX (3)  .LT.  0)  JX(3)  =  JX(3)  +  2147483399 
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c 

c  set  j x ( 1 )  =  ( ( j x (3)  +  2147483562  -  jx(2))  mod  2147483562)  +  1 
K  =  JX(3)  -  JX(2) 

IF  (K  .LE.  0)  K  =  K  +  2147483562 
c 

c  put  the  combination  back  into  jx(l) 

JX(1)  =  K 
c 

c  put  it  on  the  interval  (0.,1.) 

EUNF  =  K*4.656613E-10 

RETURN 

END 

^3fcs4c9fc9(c3|c34c9|c9fc9fc3|c3|c3)c34c3fca|c3|c9fc9tca^3|c9|c94ca4c3|c3f:3|c9fc94c3fc3|c9|ca|c3(c9^ca|e34c9|c3|c3|e34ea|c3fe3|e3|c34c9|c3fc3fe34ca|e9fca4e^c9|c34c^c3|ca|e3|ca(ca|ca|c3fe9|c 

c 

c  THIS  IS  PROGRAM  FOR  MSP  USING  BERNOULLI  POPULATIONS, 
c  SETUP  FOR  MAX  OF  5  POPULATIONS  AND  50  VECTORS, 
c  NOTE  THIS  CODE  BREAKS  ALL  TIES  WITHIN  A  VECTOR  BY  A  RANDOM  DRAW 
c  RATHER  THAN  ASSIGNING  A  FRACTION  OF  A  WIN. 
c 

c  BERNOULLI  PARAMETERS:  PX  =  BERNOULLI  SUCCESS  PROBABILITY  FOR 
c  FIRST  POPULATION  (ASSUMED  TO  BE  BEST) ,  PY  =  BERNOULLI  SUCCESS 
c  PROBABILITY  FOR  ALL  REMAINING  POPULATIONS, 
c 

PROGRAM  MNBERN 
c 
c 

REAL  PCS1(50) ,PCS2(50) ,SE1(50) ,SE2(50) 

PARAMETER (PX= . 5 ,PY= . 4264 ,NREPS=10000) 

INTEGERS  JX1(3) , JX2(3) , JX3(3) , JX4(3) , JX5(3) ,LX(3) ,MX(3) 

REAL  X(5,50) ,P(5) ,NWIN1(5,50) ,NWIN2(5,50) 
c 

c  FILES  TO  OUTPUT  RESULTS  TO 
c 

OPEN (4 , FILE= ’ bk5t 1 2 . dat  \  STATUS= » UNKNOWN ’ ) 
c 

NP0P=5 
NVECT=50 
DO  10  1=1,50 
PCS1(I)=0 
PCS2 (I)=0 
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10  CONTINUE 
c 

c  INITIALIZE  VARIABLES 
c 

c  SEEDS  FOR  L’ECUYER’S  RN  GENERATOR 
c 

JX1 (2)=748932582 
JX1 (3)=639287 
JX2 (2) =64298628 
JX2(3)=196998 
JX3 (2) =40689408 
JX3 (3) =122595 154 
JX4 (2) =680620100 
JX4 (3) =3350831 18 
JX5 (2) =64918046 
JX5 (3) =812715188 
LX (2) =427126672 
LX (3) =699944973 
MX (2) =257567734 
MX (3) =189120895 
DO  2000  1=1 , NREPS 
c 

DO  30  K=1 , NVECT 
DO  20  J=1,NP0P 
NWIN1 ( J ,K)=0 
NWIN2 ( J ,K)=0 
20  CONTINUE 

30  CONTINUE 

c 

c  RANDOM  VARIATE  GENERATION  DONE  HERE, 
c 

82  DO  87  K=1 ,50 

IF (RUNF( JX1) . LE.PX)THEN 
X(1 ,K)=1 .0 
ELSE 

X(1 ,K)=0 . 0 
END  IF 

IF(RUNF( JX2) . LE . PY) THEN 
X(2,K)=1 .0 
ELSE 

X(2,K)=0.0 
END  IF 
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IF(NP0P.EQ.2)G0  TO  87 
IF(RUNF(JX3) . LE.PY)THEN 
X(3,K)=1 .0 
ELSE 

X(3,K)=0.0 
END  IF 

IF(NP0P.EQ.3)G0  TO  87 
IF (RUNF ( JX4) . LE.PY)THEN 
X(4,K)=1 .0 
ELSE 

X(4,K)=0.0 
END  IF 

IF(NP0P.EQ.4)G0  TO  87 
IF(RUNF(JX5) .LE.PY)THEN 
X(5,K)=1 .0 
ELSE 

X(5,K)=0.0 
END  IF 

87  CONTINUE 
c 

c  METHOD  I  (BEM) 

c  SIMPLY  NEED  TO  COUNT  NUMBER  OF  TIMES  EACH  Xi  IS  BEST  IN 

c  SAME  VECTOR.  NOTE  WITH  USE  OF  BERNOULLI  POPULATIONS  HAVE 

c  INTRODUCED  ALOT  OF  TIES  THAT  MUST  BE  DEALT  WITH, 

c 

c  THIS  SECTION  ADDED  TO  TAKE  CARE  OF  TIES  WITH  BERNOULLI  POPS, 
c  SUM  IS  ACROSS  EACH  VECTOR.  BASIC  LOGIC  IN  THIS  SECTION  IS  IF 
c  THERE  ARE  ANY  TIES,  DRAW  A  RANDOM  NUMBER  (PRN  BELOW)  AND  ASSIGN 
c  FIRST  FRACTION  (0  -  1/SUM)  BASED  ON  NUMBER  OF  TIES  TO  FIRST  POP 
c  TIED  AS  WINNING  (COULD  BE  ALL  POPULATIONS  TIED  AT  ZERO)  AND 
c  NEXT  FRACTION  (1/SUM-2/SUM)  TO  NEXT  WINNING  POPULATION,  ETC. 
c 

c  10  -  FLAG  SET  TO  1  IF  SUM  IS  ZERO,  LEFT  AT  ZERO  OTHERWISE 
c  N1  -  COUNTER  FOR  NUMBER  OF  POPULATIONS  WITH  A  »1»  IN  A  VECTOR 
c  PRN  -  RANDOM  NUMBER  USED  TO  PICK  SINGLE  BEST  POP  FOR  VECTOR 
c 

101  DO  106  K=1 ,NVECT 
SUM=0 
10=0 
N1=0 

DO  103  L=1 ,NPOP 
SUM=SUM+X(L,K) 
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103 


104 

106 


107 

108 


CONTINUE 
IF (SUM . EQ . 0)THEN 
SUM=RE AL ( NP  0  P ) 

10=1 
END  IF 
PRN=RUNF(LX) 

IF ( SUM . EQ . 1 ) PRN=0 
DO  104  L=1 ,NPOP 

IF(X(L ,K) . EQ . 0 . AND . 10 . EQ . 0)G0  TO  104 
N1=N1+1 

IF(PRN.LT. (N1/SUM))THEN 
NWIN1(L,K)=NWIN1 (L,K)+1 
GO  TO  106 
END  IF 
CONTINUE 
CONTINUE 
DO  108  K=2 ,NVECT 
DO  107  J=1 ,NPOP 

NWIN1 ( J ,K) =NWIN1 ( J ,K)+NWIN1 ( J ,K-1) 
CONTINUE 
CONTINUE 


c 

c  IF  BEST  POPULATION  (#1)  WINS  OUTRIGHT,  INCREASE  PCS  BY  ONE. 
c  IF  ANOTHER  POPULATION  WINS,  DO  NOT  INCREASE  PCS. 

C  FOR  TIES  INCLUDING  BEST  POPULATION,  NEED  TO  KEEP  TRACK  OF  THE 

c  NUMBER  OF  TIES  AND  INCREASE  PCS  BY  1/NTIES. 


c 

DO  111  K=1 ,NVECT 
NTIES=1 

DO  110  J=2 ,NPOP 

IF(NWIN1(1,K) .LT.NWIN1(J,K))G0  TO  111 
IF (NWIN1 ( 1 , K) . EQ . NWIN1 ( J , K) ) NTIES=NTIES+1 
IF ( J . EQ . NPOP) PCS1 (K)=PCS1 (K) +1 . O/REAL (NTIES) 

110  CONTINUE 

111  CONTINUE 
c 

c  METHOD  II  (AVC) 

c  COMPARE  EACH  Xi  WITH  EVERY  POSSIBLE  VECTOR  V**K 

c  COMPARISONS.  THE  LOOP  INDEXES  GET  RATHER  CONFUSING 

c  HERE,  BUT  ACTUALLY  FORM  ALL  V**K  POSSIBLE  PSEDUO-VECTORS 

c  FOR  AVC.  TO  PARALLEL  BEM  IN  TERMS  OF  BREAKING  TIES  NEED 

c  TO  CONSTRUCT  EACH  PSEUDO-VECTOR. 
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c 


GO  T0(200,300,400,500) (NPOP-1) 


c  LOOPS  FOR  NP0P=2 
c 

200  DO  291  NS=1,NVECT 
DO  290  L2=1,NS 
DO  220  Ll=l ,NS 

IF(L2.LT. NS. AND. LI. LT.NS)GO  TO  220 
SUM=0 
10=0 
N1=0 

P ( 1 ) =X ( 1 ,L2) 

P(2)=X(2,L1) 

SUM=P(1)+P(2) 

IF ( SUM . EQ , 0 ) THEN 
SUM=REAL (NPOP) 

10=1 
END  IF 
PRN=RUNF (MX) 

IFCSUM.EQ . 1)PRN=0 
DO  210  L=1 ,NPOP 

IF (P (L) . EQ . 0 . AND . 10 . EQ . 0) GO  TO  210 
N1=N1+1 

IFCPRN.LT. (N1/SUM))THEN 

NWIN2 (L , NS) =NWIN2 (L , NS) +1 
GO  TO  220 
END  IF 

210  CONTINUE 

220  CONTINUE 

290  CONTINUE 

291  CONTINUE 

GO  TO  1100 
c 

c  LOOPS  FOR  NP0P=3 
c 

300  DO  391  NS=1 ,NVECT 
DO  390  L3=l ,NS 
DO  330  L2=l ,NS 
DO  320  Ll=l ,NS 

IFCL3.LT. NS. AND.L2.LT. NS. AND. LI. LT.NS)GO  TO  320 
SUM=0 
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10=0 

N1=0 

P(1)=X(1 ,L3) 

P(2)=X(2,L2) 

P(3)=X(3,L1) 

SUM=P(1)+P(2)+P(3) 

IF(SUM.EQ.0)THEN 
SUM=REAL (NPOP) 

10=1 
END  IF 
PRN=RUNF (MX) 

IF (SUM.EQ . 1)PRN=0 
DO  310  L=1 ,NPOP 

IF(P(L) .EQ.O.AND. IO.EQ.O)GO  TO  310 
N1=N1+1 

IF (PRN . LT . (Nl/SUM) ) THEN 

NWIN2(L,NS)=NWIN2(L,NS)+1 
GO  TO  320 
END  IF 

310  CONTINUE 

320  CONTINUE 
330  CONTINUE 

390  CONTINUE 

391  CONTINUE 

GO  TO  1100 
c 

c  LOOPS  FOR  NP0P=4 
c 

400  DO  491  NS=1 ,NVECT 
DO  490  L4=l ,NS 
DO  440  L3=l ,NS 
DO  430  L2=l ,NS 
DO  420  L1=1,NS 

IF(L4.LT. NS. AND.L3.LT. NS. AND. 

+  L2 . LT . NS . AND . LI . LT . NS) GO  TO  420 

SUM=0 
10=0 
N1=0 

P(1)=X(1 ,L4) 

P(2)=X(2,L3) 

P(3)=X(3,L2) 

P(4)=X(4,L1) 
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SUM=P(1)+P(2)+P(3)+P(4) 

IF(SUM.EQ.O)THEN 

SUM=REAL(NPOP) 

10=1 
END  IF 
PRN=RUNF (MX) 

IF (SUM .EQ . 1)PRN=0 
DO  410  L=1 ,NPOP 

IF(P(L) .EQ.O.AND. IO.EQ.O)GO  TO  410 
N1=N1+1 

IF(PRN.LT. (Nl/SUM) )THEN 

NWIN2(L,NS)=NWIN2(L,NS)+1 
GO  TO  420 
END  IF 


410 

CONTINUE 

420 

CONTINUE 

430 

CONTINUE 

440 

CONTINUE 

490 

CONTINUE 

491 

CONTINUE 

GO  TO 

1100 

c  LOOPS  FOR  NP0P=5 
c 

500  DO  591  NS=1 ,NVECT 
DO  590  L5=l ,NS 
DO  580  L4=l ,NS 
DO  570  L3=l ,NS 
DO  530  L2=l ,NS 
DO  520  Ll=l ,NS 

IF(L5. LT.NS.AND.L4.LT. NS. AND.L3.LT. NS. 
+  AND.L2.LT. NS. AND. LI. LT.NS)GO  TO  520 

SUM=0 
10=0 
N1=0 

P(1)=X(1 ,L5) 

P(2)=X(2 ,L4) 

P(3)=X(3 ,L3) 

P(4)=X(4,L2) 

P(5)=X(5,L1) 

SUM=P(1)+P(2)+P(3)+P (4)+P(5) 
IF(SUM.EQ.O)THEN 
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SUM=REAL (NPOP) 

10=1 
END  IF 
PRN=RUNF (MX) 

IF(SUM.EQ.1)PRN=0 
DO  510  L=1 ,NPOP 

IF(P(L) .EQ.O.AND.IO.EQ.O)GO  TO  510 
N1=N1+1 

IF(PRN.LT. (N1/SUM))THEN 

NWIN2(L,NS)=NWIN2(L,NS)+1 
GO  TO  520 
END  IF 

510  CONTINUE 

520  CONTINUE 

530  CONTINUE 
570  CONTINUE 
580  CONTINUE 

590  CONTINUE 

591  CONTINUE 

1100  CONTINUE 

DO  1105  K=2 , NVECT 
DO  1102  J=1 ,NPOP 

NWIN2 ( J ,K)=NWIN2 ( J ,K)+NWIN2 ( J ,K-1) 

1102  CONTINUE 

1105  CONTINUE 

DO  1111  K=l, NVECT 
NTIES=1 

DO  1110  J=2 ,NPOP 

IF(NWIN2(1 ,K) .LT.NWIN2(J,K))G0  TO  1111 
IF (NWIN2 ( 1 , K) . EQ . NWIN2 ( J , K) ) NTIES=NTIES+1 
IF( J .EQ .NP0P)PCS2(K)=PCS2(K)+1 . O/NTIES 

1110  CONTINUE 

1111  CONTINUE 
2000  CONTINUE 

DO  2015  1=1, NVECT 

PCS2(I)=PCS2(I)/REAL(NREPS) 

PCS1 (I)=PCS1 (I)/REAL(NREPS) 

SE1 (I)=(PCS1 (I)*(l . 0-PCS1 (I) ) /REAL(NREPS) )  ** .  5 
SE2(I)=(PCS2(I)*(1 . 0-PCS2(I))/REAL(NREPS))** .5 
RT0=PCS2(I)/PCS1(I) 
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c  FILE  4  SETUP  TO  READ  DIRECTLY  INTO  MINITAB  OR  OTHER  SPREADSHEET 
c 

WRITE (4 , 2050) NPOP , I ,PCS1 (I) ,PCS2 (I) ,RT0,SE1(I) ,SE2(I) 
c 

2015  CONTINUE 

2050  FORMAT (IX , 214 , 2X , 5 (F8 . 6 , 2X) ) 

STOP 

END 

c 

^3^34c3|c3fc3|c3|c3feafe3fc3|c3|c3fc3|c3|c3fc3fc9|c3|c3|c34e3|c3tc3|c3fc9fc?|e9fc94cafc3|c3|c3|e34c3fc3|e3|e9|cafc3|c3fc3|c3|e3|e3|ea|ca|c3|ca4c3tc9|c3|c9|e34c94ca(e3|e34c3|c9tca|c^e3fc9|es|c 

(^^^^^3fc3|e3(c3|c3|e3te^3(c3|e3|ca(c3|c3fc3|c3f:a|c3|ca4e3|c4c34e?|c34c3f;3|c3|c3|c9|c3|c3(c3te9fea|c3|ca|e34c3|cafe3fc9|e3|ea|ca|c3fc3fe3Kafc3|c34e9te34ca|e?|c3|ca|c3|e3tca|c3|c9|c 

c 

c  PROGRAM  TO  GENERATE  STANDARD  MLE  ESTIMATES  FOR  MULTINOMIAL 
c  SUCCESS  PROBABILITIES  AND  NEW  AVC  ESTIMATES.  GATHERING 
c  DATA  TO  COMPARE  THE  VARIANCES  BETWEEN  THE  TWO  ESTIMATES 
c  AND  TO  CONSTRUCT  CONFIDENCE  INTERVALS, 
c  SETUP  FOR  K=2 ,3 ,  AND  5  POPULATIONS  UP  TO  V=200. 
c  INITIALLY  SET  UP  TO  ESTIMATE  ONLY  PI  (BEST  POP) , 
c  ALTHOUGH  VARIABLES  ARE  DIMENSIONED  TO  CALCULATE  ALL  Pj’S. 
c  ONLY  MODELS  EXPONENTIAL  POPULATIONS, 
c 

c  VARIABLE  TABLE 

c 

c  NOTE:  ONLY  VARIABLES  UNIQUE  TO  THIS  PROGRAM  DEFINED  HERE, 
c  SEE  VARIABLE  TABLE  IN  PROGRAM  MSP  FOR  MORE, 

c 

C  NAME 

c 

c  AP ( J) 

c 

c  BP ( J) 

c 

c  VK 

c 

c  PA( J) 

c 

c  NLl(I) 

c 
c 


DESCRIPTION 

AVC  ESTIMATE  FOR  Pj  (P-BARj ) . 

BEM  (OR  MLE)  ESTIMATE  FOR  Pj  (P-HATj ) . 
V**K . 

ACTUAL  Pj  FOR  POPULATION  J. 

COUNT  USED  IN  CALCULATING  FIRST  TYPE  OF 
U-STAT  COVARIANCE  TERM  -  L(1,0).  INDEXED 
BY  V. 
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c 

c  NL1T  SUM  OF  NL1(I),  1=1,..., V. 

c 

c  NL2(I)  COUNT  USED  IN  CALCULATING  SECOND  TYPE  OF 

c  U-STAT  COVARIANCE  TERM  -  L(0,1). 

c 

c  NL2T  SUM  OF  NL2(I) ,  1=1, . . . ,V. 

c 

c  NV  TOTAL  NUMBER  OF  PSEUDO-VECTOR  PAIRS  WITH  A 

c  SINGLE  COMMON  ELEMENT  FOR  GIVEN  K  AND  V. 

c  USED  IN  U-STAT  COVARIANCE  CALCUALTIONS . 

c 

c  VAP(J)  VARIANCE  OF  P-BARj . 

c 

c  VBP(J)  VARIANCE  OF  P-HATj . 

c 

c  AHW  HALF-WIDTH  OF  AVC  CONFIDENCE  INTERVAL, 

c 

c  BHW  HALF-WIDTH  OF  MLE  CONFIDENCE  INTERVAL, 

c 

c  ALL(J)  AVC  LOWER-CONFIDENCE  LIMIT  FOR  P-BARj . 

c 

c  AUL(J)  AVC  UPPER-CONFIDENCE  LIMIT  FOR  P-BARj . 

c 

c  BLL(J)  MLE  LOWER-CONFIDENCE  LIMIT  FOR  P-HATj . 

c 

c  BUL(J)  MLE  UPPER-CONFIDENCE  LIMIT  FOR  P-HATj . 

c 

c  AHIT(J)  NUMBER/PCT  OF  AVC  C.I.  THAT  CAPTURE  PA(J). 

c 

c  BHIT(J)  NUMBER/PCT  OF  MLE  C.I.  THAT  CAPTURE  PA(J) . 

c 

c  WA( J)  AVERAGE  AVC  C.I.  WIDTH  FOR  P-BARj. 

c 

c  WB( J)  AVERAGE  MLE  C.I.  WIDTH  FOR  P-HATj. 

c 

c 

PROGRAM  CIK5 
c 

REAL  LAMBDA, MU, X(5, 500) ,XCNT(5) ,FC(5) 

REAL  ALL(5) ,AUL(5) ,BLL(5) ,BUL(5) ,PA(5) ,AHIT(5) ,BHIT(5) 
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REAL  XJ(5,5),NL1(200) ,NL1T,NL2(200) ,NL2T,NV,VK,NCP,NCJ 
REAL  BHW(5) , AHW(5) ,WB(5) ,WA(5) ,M3(200,200) ,M4(200,200) 
REAL  VAP(5) ,VBP(5) ,AP(5) ,BP(5) ,TAP(5) ,M5(200,200) 

REAL  NM3 , NM5 , M3MIN , M4MIN , M5MIN , M3MAX , M4MAX , M5MAX 
INTEGER*4  JX1(3) , JX2(3) , JX3(3) , JX4(3) , JX5(3) 

INTEGER  NWIN1(5) 
c 

c  FILES  TO  OUTPUT  RESULTS  TO 


OPEN (4 , FILE= ’ cik5v20a . dat ’ , STATUS= ’ UNKNOWN ’ ) 
c 

c  INITIALIZE  VARIABLES 
c 

NREPS=10000 
NVECT=20 
NP0P=5 
LAMBDA* 1.0 
MU=1 . 1227 
c 

c  ADDITIONAL  VARIABLES  TO  INITIALIZE  FOR  CONFIDENCE  INTERVAL  RUNS 
c 

Z95=l . 96 
PA(1)=.2308 
PA(2)=. 1923 
PA(3)=. 1923 
PA(4)= . 1923 
PA(5)= . 1923 
c 

c  SEEDS  FOR  L’ECUYER’S  RN  GENERATOR 
c 

JX1 (2)=748932582 
JX1 (3)=639287 
JX2 (2) =64298628 
JX2(3)=196998 
JX3 (2) =40689408 
JX3 (3) =122595154 
JX4 (2) =680620100 
JX4 (3) =3350831 18 
JX5 (2) =64918046 
JX5 (3) =812715188 
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CVL1=0 

CVL2=0 

CBP=0 

CAP=0 

CVBP=0 

CVAP=0 

BSSQ=0 

ASSQ=0 

DO  7  J=l,NPOP 
BHIT ( J)=0 
AHIT(J)=0 
WB(J)=0 
WA( J)=0 
7  CONTINUE 
10  CONTINUE 

DO  2000  1=1 , NREPS 
DO  30  J=1 ,NPOP 
BP(J)=0 
AP(J)=0 
FC( J)=0 
NWIN1 ( J)=0 
30  CONTINUE 

DO  40  K=1 ,NVECT 
NL1 (K)=l 
NL2(K)=1 

40  CONTINUE 

c 

c  RANDOM  VARIATE  GENERATION  DONE  HERE.  SEPARATE  SET  OF  RN 
c  STREAMS  FOR  EACH  POPULTION. 
c 

DO  60  K=1 ,NVECT 

X ( 1 , K) =- ( 1/LAMBDA) *LOG ( 1-RUNF ( JX1 ) ) 

X(2 ,K)=-(1/MU) *LOG ( 1-RUNF (JX2) ) 

IF(NP0P.EQ.2)G0  TO  60 
X (3 , K) =- (1/MU) *LOG ( 1-RUNF (JX3) ) 

IF(NPOP.EQ .3) GO  TO  60 
X (4 , K) =- (1/MU) *LOG (1-RUNF ( JX4) ) 

IF(NP0P.EQ.4)G0  TO  60 
X(5,K)=-(1/MU)*L0G(1-RUNF( JX5) ) 

60  CONTINUE 

c 

c  METHOD  I  (BEM  OR  MLE) 
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c  FOR  POINT  ESTIMATES  AND  CONFIDENCE  INTERVALS 

c  DO  NOT  CALCULATE  PCS. 

c 

89  DO  100  K=1 ,NVECT 

BEST=X(1 ,K) 

DO  90  J=2 ,NPOP 

TM=AMAX1(BEST,X(J,K)) 

IF (TM . GT . BEST) BEST=TM 

90  CONTINUE 

DO  95  L=1 ,NPOP 

IF(X(L,K) .EQ.TM)NWIN1(L)=NWIN1(L)+1 
95  CONTINUE 

100  CONTINUE 

DO  115  J=1 , NPOP 

BP ( J) =NWIN1 ( J) /REAL (NVECT) 

115  CONTINUE 

BPSQ=BP(1)*BP(1) 

BSSQ=BSSQ+BPSQ 

c 

c  METHOD  II  (AVC) 

c 

NM2=NVECT 

NL1T=0 

NL2T=0 

DO  900  Kl=l ,NM2 
DO  800  K2=l ,NM2 
M3(K1 ,K2)=0 
M4(K1 ,K2)=0 
M5(K1 ,K2)=0 
800  CONTINUE 

900  CONTINUE 

IF(NM2.EQ.1)THEN 

NV=1 

ELSE 

RNM2=REAL(NM2) 

RNPOP=REAL(NPOP) 

NV=RNM2**RNP0P* (RNM2-1 . 0) ** (RNPOP-1) /2 . 0 
END  IF 

VK=RNM2**RNP0P 

CALL  SORTV (NPOP ,NM2 , X) 

DO  1160  K=1 ,NM2 
DO  1148  J=1 ,  NPOP 
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XCNT( J)=l 
DO  1146  M=1 ,NPOP 
DO  1142  L=1 ,  NM2 

IF  (M . NE . J . AND . X ( J , K) . GT.X(M,L))THEN 
X J ( J , M) =NM2- (L- 1 ) 

GO  TO  1145 
ELSE 

IF(L.EQ.NM2)XJ(J,M)=0 
END  IF 

1142  CONTINUE 

1145  IF (M . NE . J) XCNT ( J ) =XCNT (J)*XJ(J,M) 

1146  CONTINUE 

1148  CONTINUE 

DO  1149  J=1 ,NPOP 

FC( J)=FC( J)+XCNT( J) 

1149  CONTINUE 
c 

c  THE  COVARIANCE  CALCULATIONS  ARE  DONE  HERE.  QUITE  DIFFERENT 
c  FOR  DIFFERENT  K,  SO  HAVE  SEPARATE  SECTIONS  BASED  ON  K.  BASIC 
c  IDEA  IS  TO  COUNT  UP  NUMBER  OF  VECTOR  PAIRS  THAT  MEET  THE 
c  CONDITIONS  FOR  THE  PROBABILITY  STATMENT  ASSOCIATED  WITH 
c  L(1 ,0)  OR  L (0 , 1) .  DO  NOT  HAVE  TO  FORM  ALL  V**K  PSUEDO-REPS 
c  TO  SOLVE.  CAN  COUNT  UP  NUMBER  OF  PAIRS  WHERE  CONDITIONS  ARE 
C  MET  AND  THEN  SUBTRACT  OFF  ANY  THAT  HAVE  COMMON  ELEMENTS  OTHER 
c  THAN  THE  ONE  DESIRED, 
c 

IF (NPOP . EQ . 2) THEN 

NLl(l)=XJ(l,2)*(XJ(l,2)-l)/2 
NLl(2)=(NM2-XJ(l,2))*(NM2-XJ(l,2)-l)/2 
IF (XJ (1 ,2) . LE . 1 ) NL 1 ( 1 ) =0 
IF (XJ (1 , 2) . EQ . NM2)NL1 (2) =0 
NL1T=NL1T+NL1 (1) 

NL1T=NL1T+NL1 (2) 

NL2(1)=(NM2-XJ (2,1))* (NM2-XJ (2,l)-l)/2 
NL2(2)=XJ(2,l)*(XJ(2,l)-l)/2 
IF(XJ (2 , 1) . EQ .NM2)NL2(1)=0 
IF (XJ (2 , 1) .LE. 1)NL2(2)=0 
NL2T=NL2T+NL2 ( 1 ) 

NL2T=NL2T+NL2 (2) 

GO  TO  1160 

END  IF 

IF (NPOP. EQ. 3) THEN 
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NC J=MAX (X J ( 1 , 2) ,XJ(1,3)) 

NCP=MIN(XJ(1,2) ,XJ(1,3)> 

NL1 ( 1) = (NCP* (NCP-1) /2) * (NCJ**2-NCJ) 

IF(XJ(1,2) . LE . 1 . 0R.XJ(1 ,3) . LE . 1)NL1 (1)=0 
NL1T=NL1T+NL1 (1) 

END  IF 

DO  1158  K2=l ,NM2 

IF (X(l ,K) . GT.X(2,K2) )M3(K,K2)=XJ(1 ,3) 

IF (NPOP .EQ . 3) GO  TO  1158 

IF(X(1 ,K) . GT.X(2 ,K2) )M4(K ,K2)=XJ(1 ,4) 

IF(X(1 ,K) . GT.X(2,K2))M5(K,K2)=XJ(1,5) 

1158  CONTINUE 

IF (NPOP. EQ. 5) THEN 
BU=XJ(1 ,2) 

BW=XJ(1 ,3) 

BY=XJ(1 ,4) 

BZ=XJ(1,5) 

WYMIN=MIN (BW , BY) 

WZMIN=MIN (BW , BZ ) 

YZMIN=MIN (BY,BZ) 

NL1  ( 1 )  =BU*  (BU-1)  /2*  ( (BW*BY*BZ)*=t=2- (BW*BY+BW*BZ+BY*BZ 
+  -WYMIN-WZMIN-YZMIN) *BW*BY*BZ) 

NL1T=NL1T+NL1 ( 1 ) 

END  IF 

1160  CONTINUE 

IF (NPOP. EQ. 2) GO  TO  1168 

NM3=0 

NM5=0 

DO  1165  Kl=l ,NM2 
DO  1163  K3=l ,NM2 

IF(K3.LE.K1)G0  TO  1163 
DO  1162  K2=l ,NM2 

IF (NPOP .EQ . 5)G0  TO  1161 

NM3=M3 (K1 , K2 ) *M3 (K3 , K2 ) -MIN (M3 (K1 , K2 ) , M3 (K3 , K2) ) 

NL2T=NL2T+NM3 

GO  TO  1162 

1161  M51=M3(K1 ,K2)*M4(K1 ,K2)*M5(K1 ,K2) 

M52=M3 (K3 , K2) *M4 (K3 , K2) *M5 (K3 , K2) 
M3MIN=MIN(M3(K1 ,K2) ,M3(K3,K2)) 

M4MIN=MIN(M4(K1 ,K2) ,M4(K3,K2)) 

M5MIN=MIN(M5(K1 ,K2) ,M5(K3,K2)) 

M3MAX=MAX(M3(K1 ,K2) ,M3(K3,K2)) 
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M4MAX=MAX(M4(K1 ,K2) ,M4(K3 ,K2) ) 

M5MAX=MAX(M5(K1,K2) ,M5(K3,K2)) 
A1=M3MAX*M4MAX-MIN(M3MAX,M4MAX)+1 
A2=M3MAX*M5MAX-MIN (M3MAX , M5MAX) 

A3=M4MAX*M5MAX-MIN (M4MAX , M5MAX) 

NM5=M51*M52- (A1+A2+A3) *M3MIN*M4MIN*M5MIN 
NL2T=NL2T+NM5 

1162  CONTINUE 

1163  CONTINUE 
1165  CONTINUE 

1168  CONTINUE 

DO  1169  J=1 ,NPOP 
TAP( J)=FC(J)/VK 
AP(J)=AP(J)+TAP(J) 

1169  CONTINUE 
APSQ=AP(1)*AP(1) 

ASSQ=ASSQ+APSQ 

CVL1=NL1T/NV 

CVL2=NL2T/NV 

c 

c  THE  FOLLOWING  SECTION  DOES  FINAL  CALCULATIONS  FOR  Pj’S  AND  Cl’S 
c 

CVL1=CVL1-AP ( 1 ) **2 
CVL2=CVL2-AP ( 1 ) **2 
VAP ( 1 ) = (CVL1+ (NPOP- 1 ) *CVL2) /NVECT 
VBP ( 1 )  =  (BP ( 1 ) *  ( 1 . O-BP ( 1 ) ) ) /NVECT 
c 

c  THESE  ARE  JUST  RUNNING  TOTALS  FOR  POINT  ESTIMATES  AND 
c  VARIANCES 
c 

CBP=CBP+BP(1) 

CAP=CAP+AP(i) 

CVBP=CVBP+VBP ( 1 ) 

CVAP=CVAP+VAP (1) 
c 

c  HERE  IS  WHERE  THE  FINAL  CONFIDENCE  INTERVAL  CONSTRUCTION  IS 
c  DONE.  THE  ONLY  DATA  WE  ARE  GATHERING  HERE  IS  WHETHER  OR  NOT 

c  EACH  AVC  AND  BEM  Cl  IS  A  HIT  OR  NOT,  AND  THE  WIDTH  OF  THE 

c  INTERVAL . 
c 

DO  1200  J=l,l 

BHW ( J) =Z95*SQRT (VBP ( J) ) 
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BLL ( J ) =BP ( J) -BHW ( J ) 

BUL( J)=BP(J)+BHW(J) 

WB( J)=WB( J)+2*BHW( J) 

IF (BLL( J) .LE.PA(J) . AND . BUL ( J ) .GE.PA(J) )BHIT( J)=BHIT(J)+1 
AHW(J)=Z95*SQRT(VAP(J) ) 

ALL( J)=AP( J) -AHW( J) 

AUL(J)=AP(J)+AHW(J) 

WA( J)=WA( J)+2*AHW( J) 

IF(ALL(J) .LE.PA(J) . AND. AUL(J) .GE.PA( J) ) AHIT(J)=AHIT(J)+1 
1200  CONTINUE 
2000  CONTINUE 
c 

c  HERE  WE  CALCULATE  THE  PERCENTAGE  OF  HITS  WE  HAD  WITH  EACH  METHOD 
c  OF  CONFIDENCE  INTERVALS  AND  ALSO  THE  AVERAGE  Cl  WIDTH 
c 

DO  2200  J=l,l 

BHIT ( J) =BHIT ( J) /REAL (NREPS) 

WB ( J ) =WB ( J ) /RE AL ( NREPS ) 

AHIT ( J) =AHIT ( J) /REAL (NREPS) 

WA( J)=WA( J) /REAL ( NREPS ) 

2200  CONTINUE 
c 

c  AS  A  CHECK  ON  THE  BIAS  OF  OUR  AVC  VARIANCE  ESTIMATE,  CALCULATE 
C  THE  SAMPLE  VARIANCE  ASSOCIATED  WITH  EACH  METHOD  HERE, 
c 

SVARB= (BSSQ- (CBP*CBP/NREPS) ) / (NREPS-1) 

SVARA= (ASSQ- (CAP*CAP/NREPS) ) / (NREPS-1) 

CBP=CBP/NREPS 

CAP=CAP/NREPS 

CVBP=CVBP/NREPS 

CVAP=CVAP/NREPS 

VBIASB=CVMP-SVARB 

VBIASA=CVAP-SVARA 

c 

c  FILE  4  SETUP  TO  PRINT  RESULTS  IN  READABLE  FASHION 
c 

WRITE(4,*) ’CONFIDENCE  INTERVAL  RESULTS  EXP  POPULATIONS’ 

WRITE (4,*) ’K=5  AND  Pl=0.2308’ 

WRITE (4,*) 

WRITE (4,*) ’BEM  RESULTS  ’ 

WRITE (4,*) ’  V  AVG  VAR  SAMPLE  VAR  VAR  BIAS’ 

WRITE (4 , 2500) NVECT , CVBP , SVARB , VBIASB 
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WRITE (4,*)’  ’ 

WRITE (4,*) ’  V  PCT  HITS  AVG  Cl  WIDTH’ 

WRITE (4 , 2600) NVECT , BHIT ( 1 ) , WB ( 1 ) 

WRITE (4,*) ’  ’ 

WRITE (4,*) ’AVC  RESULTS’ 

WRITE (4,*) ’  V  AVG  VAR  SAMPLE  VAR  VAR  BIAS’ 
WRITE (4 , 2500) NVECT , CVAP , SVARA , VBIASA 
WRITE (4,*) ’  ’ 

WRITE (4,*)’  V  PCT  HITS  AVG  Cl  WIDTH’ 

WRITE (4 , 2600) NVECT , AHIT ( 1 ) , WA ( 1 ) 

2500  FORMAT ( IX , 14 , 2X , 3 (F8 . 4 , 4X) ) 

2600  FORMAT ( IX , 14 , 2X , 2 (F8 . 4 , 4X) ) 

STOP 

END 
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